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INTRODUCTION 


1.1  Project  Description 

This  is  the  First  Yearly  Technical  Report  on  a  three  year 
effort  to  study  physical  processes  of  relevance  to  the  mass 
spectrometric  measurement  of  stratospheric  ions.  The  effort 
involves  the  development  of  a  Monte  Carlo  model  of  the  free jet 
expansion  occurring  within  the  mass  spectrometer  including  the 
effects  of  agglomeration  onto,  and  fragmentation  of,  ionic 
clusters. 

The  attempt  to  carry  out  in  situ  mass  spectrometry  in  the 
stratosphere  is  complicated  by  changes  that  may  occur  in  the 
gas  stream  as  it  expands  after  passage  through  the  orifice. 

Both  positive  and  negative  ions  exist  in  the  stratosphere  with 
clustered  polar  molecules  surrounding  the  ion  core.  As  these 
ion  clusters  are  carried  along  in  the  expanding  gas  stream,  the 
falling  temperature  will  tend  to  favor  the  formation  of  larger 
clusters.  The  charge-dipole  interaction  is  characterized  by 
large  cross  sections,  so  agglomeration  of  polar  molecules  may 
change  the  cluster  size  distribution  that  the  quadrupole  sees 
from  the  distribution  that  exists  in  the  undisturbed  stratos¬ 
phere.  Conversely,  tjie  measured  cluster  size  distribution  may 
be  driven  towards  smaller  clusters  via  fragmentation.  As  the 
ionic  clusters  are  selectively  accelerated  by  the  electric  field 
within  the  mass  spectrometer,  high  energy  collisions  with 
neutrals  may  break  apart  the  clusters. 

The  present  effort  involves  a  Monte  Carlo  simulation  of 
these  processes,  so  that  a  model  can  be  used  to  relate  the 
measured  properties  to  those  existing  in  the  undisturbed 
atmosphere.  Sections  2  through  8  describe  the  method  in 


general  terms,  although  most  of  the  new  results  were  in  fact 
accomplished  in  support  of  the  present  contract.  Section  9 
describes  particular  results  obtained  when  the  procedures  were 
applied  to  the  freejet  expansion  within  a  mass  spectrometer. 

1.2  Overview  of  the  Direct  Simulation  Monte  Carlo  Method 

The  direct  simulation  Monte  Carlo  method,  as  pioneered  by 
G.  A.  Bird,^  provides  a  powerful  technique  for  the  simulation 
of  real  gas  flows.  It  bridges  the  gap  between  continuum  and 
free  molecular  flow,  retaining  validity  in  either  extreme.  It 
can  be  used  to  describe  complex  mixtures,  including  effects  of 
chemical  reactions,  heat  conduction,  viscosity  and  diffusion 
for  flows  in  three  dimensions.  To  date,  it  is  the  only 
approach  which  can  demonstrate  these  abilities  for  general 
flow  configurations. 

The  basic  calculational  technique  is  well  described  by  its 
originator  in  Reference  1,  to  which  frequent  reference  will  be 
made.  The  present  purpose  is  to  describe  how  the  technique  is 
implemented  at  Spectral  Sciences,  Inc.  (SSI),  with  special 
emphasis  on  extensions  developed  at  SSI  and  elsewhere  after 
the  publication  of  Reference  1.  Elementary  concepts  and 
relations  which  are  essential  to  a  coherent  explanation  are 
included  here  for  clarity. 

The  direct  simulation  Monte  Carlo  method  involves  storing 
a  discrete  number  of  molecules  (via  their  velocities,  positions 
and  other  pertinent  information)  in  a  computer.  The  solution 


t - 

Bird,  G.  A.,  Molecular  Gas  Dynamics,  Clarendon  Press, 
Oxford,  1976. 
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region  is  broken  up  into  a  number  of  separate  cells,  and  the 

solution  is  stepped  forward  in  time  in  a  two  stage  process. 

First,  the  molecules  are  advanced  along  their  trajectories  by 

an  amount  appropriate  to  their  velocity  and  a  time  increment 

At  .  In  this  first  stage  some  molecules  will  leave  the  solution 
m 

region  and  some  will  be  introduced  as  determined  by  the  boundary 
conditions  for  a  particular  problem.  The  second  stage  is  to 
simulate  collisions  in  each  cell  appropriate  to  Atm  so  that 
collision  frequencies  are  properly  simulated.  A  basic  hypothe¬ 
sis  of  the  method  is  that  if  the  time  step  is  made  small  enough 
the  processes  of  translations  and  collisions  can  be  uncoupled 
in  this  manner. 

Periodically,  the  solution  is  sampled  by  accumulating 
statistical  sums  of  number  densities,  velocities  and  other 
basic  properties.  The  solution  is  run  repeatedly  until  statis¬ 
tical  deviations  are  reduced  to  a  desired  limit,  and  then 
physically  meaningful  output  quantities  are  computed  from  the 
statistical  sums.  The  number  of  molecules  represented  is  typi¬ 
cally  a  few  thousand  at  a  time,  which  is  vastly  fewer  than  the 
number  of  molecules  occurring  in  virtually  all  real  flows. 

Hence,  the  construction  of  a  dynamically  similar  flow  to  be 
simulated  in  the  computer  is  an  essential  feature  of  the 
method. 

2.  GAS  MODEL  AND  EQUILIBRIUM  PROPERTIES 
2.1  Preliminary  Equilibrium  Gas  Relations 

For  most  problems  of  interest  there  is  a  far  field  equil¬ 
ibrium  state  whose  properties  are  of  relevance  to  the  problem 
to  be  solved.  Frequently  length  and  velocity  scales  for  the 
problem  are  obtained  from  the  far  field  state  and  used  to 


nondimensionalize  the  internal  code  variables.  Even  when  the 
far  field  state  is  not  used  for  scaling  purposes,  it  still 
provides  an  important  comparison  case. 


For  a  rest  gas  in  equilibrium  the  normalized  distribution 

function  for  the  relative  speed,  c  ,  between  molecules  of  spe- 

r  2 

cies  i  and  molecules  of  species  j  is  given  by 


4a?/2 


fij  «=r> 


exp 


KjCr) 
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where 

a  -  J!ii. 

ij  "  ' 

and  \i..  is  the  reduced  mass  of  the  pair;  i.e.. 


(2! 


m^nij 

yij  "  m.^  +  m^ 


(3] 


with  m^  and  m^  representing  the  masses  of  the  two  species.  In 
these  relations,  is  the  far  field  temperature  and  RQ  is  the 
universal  gas  constant.  (Rq  is  used  instead  of  Boltzmann's 
constant  since  the  molecular  masses  will  be  consistently  repre¬ 
sented  in  atomic  mass  units  rather  than  grams.)  The  available 
translational  collision  energy  between  two  molecules,  E  ,  is 
given  by 


(4 
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Chapman,  Sydney,  and  Cowling,  T.  G.,  The  Mathematical  Theory 
of  Non-Uniform  Gases,  ’rd  ed.  Cambridge  University  Press, 


2.2  Analytical  Form  of  The  Collision  Cross  Section 


Whenever  the  direct  simulation  Monte  Carlo  method  is 
applied,  it  is  necessary  to  make  tradeoffs  between  accuracy 
and  simplicity  in  molecular  models.  It  does  no  good  to  use  a 
complex  molecular  potential  surface  and  then  find  that  reason¬ 
able  computer  run  times  result  in  very  large  statistical  fluc¬ 
tuations  in  the  output.  Since  the  final  output  will  reflect 
errors  in  the  statistics  as  well  as  errors  in  the  models,  there 
is  a  strong  impetus  to  use  models  which  contain  the  essential 
physics,  but  which  can  be  applied  in  a  computationally  efficient 

manner.  The  current  state-of-the-art  is  the  Variable-Hard- 

3 

Sphere  (VHS)  model.  In  this  model  molecules  have  a  collision 
cross  section  which  varies  as  an  inverse  power  of  the  available 
collision  energy.  Hence,  if  is  the  collision  cross  section 
for  collisions  of  species  i  with  species  j,  then  ck..  is  given 
by  a  relation  of  the  form 


-u> 


ij 


(5) 


where  is  a  constant  coefficient.  It  follows  that  the 
effective  diameter  for  molecules  of  species  i,  d.,  is  implicitly 
defined  as  a  function  of  available  collision  energy  by  the 
relation 


xx 
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Aii  can  ke  determined  from  a  reference  cross  section  and 
velocity  via 


(6) 


Bird,  G.  A.,  "Monte-Carlo  Simulation  in  an  Engineering  Context, 
Proceedings  of  the  12th  International  Symposium  on  Rarefied  Gas 
Dynamics,  Vol.  74,  Progress  in  Astronautics  and  Aeronautics, 
AIAA,  New  York,  198lT"^ 
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If  a  reference  cross  section  is  given  for  a  reference  tempera¬ 
ture  rather  than  a  reference  velocity,  then  the  usual  choice  for 
the  reference  velocity  is  that  velocity  which  has  a  collision 
energy  equal  to  the  mean  collision  energy  occurring  in  colli¬ 
sions  at  the  reference  temperature.  Mathematically,  this  is 
equivalent  to 


(4 
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r  11 
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r  11 


where  the  bars  over  the  quantities  indicate  averages  taken  over 
the  distribution  function  given  in  Eq.  (1)  evaluated  for  m.  =  m. 
and  Tb  =  T  Equation  (8)  can  be  simplified  to  give 


(4, 


4(2  -  “’Vref 


For  simulations  involving  a  large  number  of  species,  ref¬ 
erence  cross  sections  are  frequently  not  available  for  all 
possible  collision  pairs.  In  this  case  it  is  possible  to 
specify  for  self  collisions  only,  and  then  use  Eq.  (6)  to 
get  a  molecular  diameter  as  a  function  of  collision  energy. 
Then,  applying  the  relation 


.  =  "  [(di +  dj)/2 


the  coefficient  in  Eq.  (5)  for  interspecie  collisions  is  given 
by 


[#(vsir + 


X, 


For  the  internal  workings  of  a  Monte  Carlo  code,  it  is 
usually  more  convenient  to  express  the  collision  cross  section 
as  a  function  of  the  relative  collision  velocity  rather  than 
the  collision  energy.  This  is  simply  achieved  via  the 
relation 

°ij  =  BijCr2<"  '  <12 

where 


Bij  * 


—to 

•«(^)  • 


The  parameter  oo  can  be  related  to  r\,  the  exponent  of 

distance  in  an  inverse  power  intermolecular  force  law  via  the 
3 

relation 

w  =  2/ (n  -  1)  .  (14) 

Hence,  hard  sphere  molecules  (for  which  n  goes  to  infinity)  are 
represented  by  w  equal  to  zero.  There  is  a  substantial  body  of 
evidence,  however,  that  the  effective  size  of  molecules  does 
indeed  decrease  with  increasing  collision  energy,  so  a  positive 
value  of  w  is  usually  a  better  choice,  w  can  be  determined 
from  molecular  beam  data,  or  from  its  macroscopic  implications. 
For  example,  if  s  is  the  temperature  exponent  for  the  coefficient 

3 

of  viscosity,  then  it  can  be  shown  that 

s  =  a)  +  1/2  ,  (15) 

so  a  measurement  of  the  temperature  dependence  of  the  viscosity 
coefficient  serves  as  an  indirect  determination  of  w. 

In  order  to  incorporate  the  model  for  internal  energy 
transfer  to  be  discussed  in  Section  4,  it  is  necessary  that  u> 
be  assumed  the  same  for  all  interactions.  This  represents  one 
of  the  major  restrictions  in  the  current  state  of  modeling. 


Although  the  sizes  of  molecules  are  allowed  to  vary  in  the 

VHS  molecular  model  in  deciding  whether  or  not  a  collision  is 

to  occur,  when  a  collision  does  occur  the  post  collision  velocity 

components  are  computed  as  if  it  were  a  hard  sphere  collision 

(see  Section  4).  This  results  in  a  substantial  computational 

simplification  and  yet  retains  good  agreement  with  the  macro- 

3 

scopic  predictions  of  the  more  exact  model.  (See  Reference  1 
for  a  discussion  of  molecular  scattering  for  general  power  law 
potentials. ) 


2.3  Equilibrium  Reference  Properties  for  a  Multi-Component  Gas 

One  advantage  of  the  VHS  model  is  that  the  molecules  have 
a  well  defined  cross  section,  so  it  is  possible  to  define  a  mean 
free  path  without  putting  limitations  on  the  minimum  deflection 
angle  that  is  considered.  As  is  the  general  case  for  multi- 
component  gases,  however,  each  component  has  its  own  mean  free 
path,  and  the  overall  mean  free  path  for  the  mixture  must  be 
defined  as  a  weighted  average  of  the  mean  free  paths  of  the 
individual  species.  The  somewhat  cumbersome  relations  required 
to  calculate  the  overall  mean  free  path  are  given  here.  It 
should  be  noted  that  the  mean  free  path  is  calculated  only  once 
for  a  given  problem,  so  the  computational  effort  required  to 
evaluate  it  is  completely  negligible. 

An  individual  molecule  of  species  i  will  suffer  collisions 
with  molecules  of  species  j  with  a  frequency  given  by 

.j  _  «  r— (16) 


nj®  °ijcr 


where  n^  is  the  number  density  of  species  j  and  o^jCr  is  the 
average  product  of  cross  section  times  relative  velocity  for 
the  two  species,  obtained  by  integrating  over  the  distribution 
function  given  in  Eq.  (1).  When  this  operation  is  performed, 
the  result  is 


where  r  denotes  the  gamma  function. 

The  total  collision  frequency  for  an  individual  molecule 
of  species  i,  is  obtained  by  summing  Eq.  (15)  over  all 
species,  i.e. 
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j=l 


and  the  mean  free  path,  X^,  for  molecules  of  species  i  is 
given  by 
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where  cT  is  the  mean  molecular  speed  for  species  i  molecules. 
The  mean  free  path  for  the  gas  mixture,  X^,  is  then  defined  as 
the  number  density  weighted  average  of  the  X^  via 
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where  n^  is  the  total  number  density: 
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A  useful  velocity  scale  is  given  by  v  .  defined  by 


where  m  is  the  reference  mean  molecular  weight,  i.e 


(23 


vg  is  the  most  probable  molecular  speed  for  molecules  of  the 
mean  molecular  weight  at  the  reference  temperature. 


2.4  Internal  Energy  Model 

The  current  state  of  modeling  for  internal  energy  effects 

in  Monte  Carlo  flow  field  simulations  is  the  phenomenological 

4 

model  of  Borgnakke  and  Larsen.  In  this  model,  transfer  of 
energy  between  internal  and  translational  modes  is  allowed, 
but  it  is  necessary  to  assume  that  each  species  has  a  fixed 
number  of  internal  degrees  of  freedom,  This  is  equivalent 

to  assuming  a  constant  specific  heat,  Cp^,  for  each  species 
which  can  be  related  to  the  number  of  internal  degrees  of 
freedom  via 


=  2 


C  .m. 
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Alternatively,  can  be  related  to  the  ratio  of  specific  heats 


for  species  i,  by  the  relation 


(25 


The  interchange  of  internal  and  translational  energy  will  be 
discussed  in  Section  4,  and  the  selection  of  initial  conditions 
will  be  considered  in  Section  7. 


^Borgnakke,  Claus,  and  Larsen,  Paul  S.,  "Statistical  Collision 
Model  for  Monte  Carlo  Simulation  of  Polyatomic  Gas  Mixture," 
Journal  of  Computational  Physics,  Vol.  18,  1975,  pp.  405-420. 


3.  INTERNAL  REPRESENTATION 
3 . 1  State  Vector 

Each  simulated  molecule  in  the  direct  simulation  Monte 
Carlo  method  is  represented  by  a  state  vector  which  comprises 
all  of  the  information  the  code  has  with  regard  to  that  partic 
ular  molecule.  The  state  vector  has: 

•  Position  element (s)  defining  the  location  of  the 
molecule  in  the  coordinate  system  being  used. 

For  axisymmetric  simulations,  this  is  a  radial 
and  an  axial  element. 

•  Three  velocity  elements.  A  molecular  collision 
is  always  considered  as  a  three  dimensional 
event,  regardless  of  the  overall  dimensionality 
of  the  problem.  For  spatially  one  dimensional 
problems  it  is  possible  to  store  only  two  pieces 
of  velocity  information  and  compute  the  required 
three  velocity  components  as  needed  for  collision 
sampling.  This  is  an  example  of  the  frequent 
tradeoff  which  must  be  made  between  storage  and 
computing  requirements. 

•  A  value  for  the  internal  energy  of  the  molecule. 

Note  that  the  basic  model  does  not  discriminate 
between  internal  modes  for  a  particular  species. 

This  can  be  done,  if  desired,  by  introducing 
separate  species  for  the  separate  modes. 

•  An  indicator  determining  the  molecular  species. 

This  indicator  in  turn  implies  all  of  the 
properties  associated  with  that  species  (molec¬ 
ular  weight,  number  of  internal  degrees  of 
freedom,  name,  etc.). 


•  An  indicator  giving  the  cell  in  which  the  mole¬ 
cule  currently  resides.  It  is  possible  to 
avoid  allocating  this  particular  storage  element, 
but  it  usually  results  in  enough  computational 
simplification  to  justify  its  use. 

3.2  Reduction  to  a  Reasonable  Number  of  Simulated  Molecules 

It  is  clearly  impossible  to  run  a  computer  simulation  with 
anywhere  near  the  same  number  of  molecules  that  exist  in  the 
actual  flow  problem.  The  adjustment  that  is  made  to  make  the 
simulation  possible  is  to  artificially  increase  the  cross 
section,  and  decrease  the  number  density,  by  a  large  factor. 

It  is  the  product  of  number  density  and  cross  section  which 
determines  the  collision  frequency  for  a  given  molecule,  and 
it  is  the  collision  frequency  which  must  be  correctly  simulated 
if  the  correspondence  between  the  real  and  simulated  flows  is 
to  be  correct.  This  is  an  essential  feature  of  the  direct 
simulation  method  which  has  not  always  been  adequately  empha¬ 
sized.  It  means  that  the  internal  scaling  factors  do  not  pro¬ 
ceed  on  a  strictly  dimensional  basis.  For  example,  the  scaling 
factor  for  cross  sections  is  not  the  square  of  the  scaling 
factor  for  lengths. 

3.3  Internal  Scales 

Many  problems  are  more  reasonably  handled  if  the  internal 
calculations  are  carried  out  with  scaled  or  dimensionless 
values.  This  avoids  possible  problems  such  as  numerical  over¬ 
flow  which  can  cause  an  execution  time  error.  Such  errors  can 
be  particularly  insidious  and  difficult  to  locate  in  a  code 
whose  very  essence  involves  the  random  combination  of  numbers. 
Furthermore,  examination  of  scaled  values  makes  the  detection 


of  erroneous  values  easier  while  debugging  codes  since  large 
exponents  are  usually  indicative  of  an  error  when  the  variables 
are  internally  scaled.  At  SSI,  at  least,  the  output  is  produced 
in  physically  meaningful  dimensional  form.  Hence,  the  scaling 
that  is  discussed  here  is  irrelevant  (or  nearly  so)  to  the 
interpretation  of  code  output;  it  is  strictly  a  matter  of  the 
internal  representation. 

The  obvious  choices  for  length  and  velocity  scales  are  X^ 

and  v  as  defined  in  Section  2,  which  are  used  to  nondimension- 
s 

alize  the  position  and  velocity  elements  of  the  state  vector. 

There  is  no  need  to  provide  further  nondimensionalization  of 

mass  beyond  representing  them  in  atomic  mass  units,  so  none  is 

2 

provided.  Hence,  the  scaling  factor  for  energy  is  just  v  , 

S 

which  is  used  to  nondimensionalize  the  internal  energy  element 
of  the  state  vector. 

Number  densities  are  scaled  with  respect  to  the  far  field 
reference  number  density,  nw,  which  leaves  only  the  cross  sec¬ 
tion  scaling  factor  to  be  determined.  This  factor  follows 
from  the  condition  of  flow  similarity,  which  requires  that  the 
probability  of  a  molecule  suffering  a  collision  in  traveling  a 
given  path  length  be  accurately  simulated.  This  dimensionless 
probability  can  be  expressed  as  the  product  of  a  cross  section 
times  a  number  density  times  a  path  length  (at  least  for  small 
enough  path  lengths),  and  it  is  required  that  this  product  be 
the  same  for  dimesional  and  scaled  representations.  This 
implies  that  the  product  of  the  scaling  factors  for  these  three 
quantities  be  unity  and,  therefore,  the  cross  section  scaling 
factor  is  1/ (n^A^) .  The  internal  scaling  factors  used  for  the 
SSI  Monte  Carlo  codes  are  summarized  in  Table  1. 


Table  1.  Scaling  Factors  Used  for  the  Internal 
Representation  of  Quantities  in  the  SSI  Direct 
Simulation  Monte  Carlo  Codes.  All  Variables 
are  Defined  in  Section  2. 


Property 

Scaling  Factor 

Length 

Velocity 

vs 

Time 

X<*/VS 

Number  Density 

n» 

Mass 

a.m.u. 

Energy 

2 

(a.m.u. )v 

s 

Cross  Section 

1/ («•*.) 

3.4  Weighting  Factors 

Weighting  factors  are  a  crucial  element  of  a  successful 
Monte  Carlo  simulation,  allowing  trace  species  to  be  described 
with  reasonable  statistics.  The  weighting  factor  is  the  number 
of  "real"  molecules  that  corresponds  to  each  "simulated"  mole¬ 
cule.  A  "simulated"  molecule  corresponds  to  one  molecule’s 
worth  of  storage  (one  state  vector)  allocated  in  the  computer, 
and  the  weighting  factor  is  its  statistical  weight.  So,  for 
example,  the  total  number  density  in  a  cell  might  be 
represented 


ncell 


naw 

“V" 


i 


t 


(26) 


where  indicates  the  number  of  simulated  molecules  of  species 
i  in  the  cell,  is  the  weighting  factor  for  the  species  in 
that  cell,  V  is  the  cell  volume  and  p  is  the  number  of  species. 
The  product  that  appears  in  Eq.  (26)  is  termed  the  number 

of  "real"  molecules  of  species  i  in  the  cell.  Note  that  n  ^ 
as  calculated  by  Eq.  (26)  is  a  scaled  value;  it  would  have  to 
be  multipled  by  n0,  as  shown  in  Table  1,  to  become  a  dimensional 
evaluation  of  the  number  density. 

The  weighting  factors  used  in  the  SSI  codes  are  dependent 
on  cell  and  species.  Hence,  flowfields  where  a  given  species 
is  much  more  dominant  in  one  portion  of  the  solution  region 
than  another  can  be  accurately  represented.  It  is  possible, 
of  course,  to  make  weighting  factors  functions  of  other  vari¬ 
ables,  such  as  velocity,  for  specialized  purposes. 

A  critical  error  that  can  occur  in  Monte  Carlo  codes  is  to 
have  the  number  of  simulated  molecules  exceed  the  dimensioned 
limit  of  the  code.  On  the  other  hand,  it  is  generally  desir¬ 
able  to  have  as  many  molecules  as  is  feasible  to  obtain  good 
statistics.  Resolution  of  these  conflicting  considerations  is 
complicated  by  lack  of  a  priori  knowledge  of  what  the  species 
number  densities  will  be  as  a  function  of  space  and  time.  The 
way  the  resolution  is  achieved  in  the  SSI  codes  is  by  a  dynamic 
adjustment  of  the  weighting  factors,  as  required.  This  keeps 
the  number  of  simulated  molecules  more  or  less  constant  while 
allowing  the  number  of  real  molecules  to  adjust  as  the  solution 
evolves.  The  introduction  of  weighting  factors,  with  the 
ability  to  adjust  them  as  the  solution  demands,  is  an  important 
feature  of  a  successful  Monte  Carlo  description. 
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COLLISION  MECHANICS 


4.1  Relations  for  Elastic  Collisions 

The  purpose  of  this  section  is  to  present  relations 
appropriate  to  the  simulation  of  a  collision  in  a  Monte  Carlo 
flowfield  code.  (The  question  of  how  molecules  are  selected 
for  collisions,  which  is  crucial  to  the  proper  simulation  of 
collision  frequency,  will  be  taken  up  in  Section  6.)  Conser¬ 
vation  of  momentum  implies  that  the  center-of-mass  velocity  of 
the  the  collision  pair  is  unchanged  by  the  collision;  and  con¬ 
servation  of  energy  then  implies  that  the  magnitude  of  the 
relative  velocity  between  the  collision  partners  is  also 
unchanged  by  the  collision.^  Since  the  collision  is  treated 
as  a  statistical  event,  all  that  remains  is  to  select  the 
direction  of  the  post -collision  relative  velocity  vector  from 
the  correct  distribution.  As  mentioned  in  Section  2.2,  colli¬ 
sions  in  the  VHS  model  are  treated  as  hard  sphere  collisions 
when  they  occur  (though  they  do  not  occur  with  the  same  velocity 
dependence  as  do  hard  sphere  collisions).  Hence,  as  far  as  the 
collision  mechanics  is  concerned,  the  model  is  a  hard  sphere 
model.  For  hard  sphere  molecules,  all  directions  for  the  post¬ 
collision  relative  velocity  vector  are  equally  likely.  This  is 
the  chief  computational  simplicity  of  the  VHS  model. 

Let  the  two  molecules  be  identified  by  subscripts  1  and  2, 
with  m  and  v  denoting  their  masses  and  velocities.  If  i  and  f 
indicate  initial  and  final  states,  then  the  relations  for  the 
collision  can  be  summarized  via: 


Vincenti,  Walter  G.,  and  Kruger,  Charles  H.,  Jr.,  Introduction 
to  Physical  Gas  Dynamics,  John  Wiley  and  Sons,  1965,  pp.  348- 


cos  (6) 

=  1-26  , 

(29) 

sin(6) 

=  yjl  -  COS2(0)  , 

(30) 

0  = 

2  it  6  , 

(31) 

*rf  = 

vr£cos(0),  sin (0)cos (<J>) ,  sin (0) sin ( <J> )J 
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(32) 

*lf  = 

-*■  m2 

vcm  +  m^  +  mj  Vrf  * 

(33) 

*2f  = 

-v  ^ 

V  —  v  ^ 

cm  m1  +  n»2  rf  ' 

(34) 

In  these  relations,  and  throughout  this  report,  $  indicates 
a  random  variable  which  is  evenly  distributed  on  the  interval 
zero  to  one.  Each  time  that  B  appears  a  different  evaluation 
of  the  random  variable  is  implied.  Note  that  the  expression 
for  the  post-collision  relative  velocity  vector  (Eq.  (32))  is 
not  coordinate  system  specific.  The  indicated  vector  components 
can  apply  to  any  locally  orthogonal  coordinate  system,  since  the 
direction  implied  is  random.  The  convenient  coordinate  system 
to  use,  of  course,  is  the  coordinate  system  used  to  define  the 
velocity  elements  of  the  state  vector.  For  axisymmetric  simu¬ 
lations  this  will  be  radial,  azimuthal  and  axial  velocity 


components. 


1 


i 


l 


c, 

V 


V 

v* 
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4.2  Simulation  of  Inelastic  Collisions 

4 

The  SSI  codes  use  the  Borgnakke  and  Larsen  phenomenolog¬ 
ical  model  for  transfer  of  energy  between  internal  and  transla¬ 
tional  modes.  In  this  model,  a  collision  is  assumed  to  be 
either  perfectly  elastic  or  perfectly  inelastic,  via  a  user 
specified  probability.  A  perfectly  inelastic  collision  is 
achieved  by  summing  the  total  pre-collision  energy  (internal 
energy  of  both  molecules  plus  the  translational  energy  of  their 
relative  motion,  Eq.  (4))  and  then  assigning  post-collision 
values  from  the  equilibrium  distribution  for  collisions  with 
that  total  amount  of  energy,  taking  into  account  the  number  of 
internal  degrees  of  freedom  in  the  two  molecules.  Note  that 
this  model  has  the  ability  to  relax  from  a  nonequilibrium  to  an 
equilibrium  state  via  an  effective  collision  number.  The  abil¬ 
ity  to  exchange  internal  energy  in  such  a  manner  comprises  a 
significant  increase  in  capability  for  Monte  Carlo  codes  beyond 
the  previous  models  where  molecules  had  no  internal  energy. 

It  is  this  capability  which  enables  the  codes  to  realistically 
predict  the  macroscopic  effects  of  polyatomic  gas  flow. 

Let  an<*  ^2  ke  number  of  internal  degrees  of  freedom 

of  the  two  molecules  in  an  inelastic  collision,  and  E  be  the 

s 

total  collision  energy  defined  by 


=  E  .  +  E, .  +  E* . 
cr  li  2i 


(35) 


where  Ec^  is  the  initial  translational  collision  energy  defined 
by  Eq.  (4),  and  E^  and  Ej^  are  the  pre-collision  internal 
energies  of  the  two  molecules.  If  £  is  defined  by 

5  =  Ecf/Es  ,  (36) 


where  Ec^  is  the  post-collision  translational  energy,  then  £  is 
selected  according  to  the  distribution 
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where 


[■<;>  -  col1' 

L  1  -  w  J 

l<cs-=-rJ 

9 

b  =  <;>  -  1  ,  (3! 

and 

<?>  =  (^  +  s21/2  .  ( 4 1 

The  sampling  of  Eq.  (37)  is  done  via  a  technique  that  is 
used  frequently  in  Monte  Carlo  flowfield  codes,  the  acceptance” 
rejection  method.  Equation  (37)  has  been  normalized  so  that 
the  maximum  value  of  the  function  is  unity,  and  the  parameter 
of  the  distribution,  £,  varies  from  zero  to  one.  The  sampling 
is  done  as  follows: 

•  Choose  a  random  value  of  £.  I.e.,  set  £  equal 

to  6,  a  random  variable. 

•  Evaluate  the  distribution  function  for  this 
value  of  £,  and  call  it  ftest» 

•  Get  a  second  random  variable,  and  check  to  see 
if  it  is  greater  or  less  than  ftest»  If  it  is 
greater  than  ftestr  go  back  to  the  first  step 
and  repeat  the  process.  If  it  is  less  than 
^test'  then  keep  the  value  of  £  obtained  in 
the  first  step. 

Note  that  the  probability  that  any  original  value  of  £  will  be 
kept  is  proportional  to  f(£).  This  is  an  extremely  general 
technique  and,  as  such,  is  very  powerful.  It  is  not  always 
efficient,  however,  and  direct  sampling  of  distributions  is 
usually  to  be  preferred  if  it  can  be  accomplished.  In  this 
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case,  if  one  molecule  is  monatomic  and  the  other  is  diatomic, 
i.e.,  if  <£>  =  1,  then  a  direct  sampling  of  the  above  distri¬ 
bution  is  achieved,  via 


K  =  B 


[1/(2-003 


(<c>  =  1) 


Once  £,  and  therefore  the  post-collision  translation 
energy,  is  determined,  the  magnitude  of  the  post-collision 
relative  velocity  is  defined  via 

vr  =  V2E  cf/p12  *  (42) 

For  inelastic  collisions,  this  relation  takes  the  place  of  Eq. 
(28)  in  the  determination  of  the  post-collision  velocity  ele¬ 
ments  of  the  state  vectors. 

The  remainder  of  the  collision  energy  must  be  divided  up 
between  the  internal  modes  of  the  two  molecules.  If  one  of  the 
molecules  is  monatomic  (i.e.,  has  zero  internal  degrees  of  free¬ 
dom)  ,  then  all  of  the  internal  energy  goes  to  the  other  by 
default.  Otherwise,  if  x  is  defined  by 

x  =  Elf/ (Es  -  Ecf)  ,  (43) 

where  Elf  is  the  post-collision  internal  energy  of  the  first 
molecule,  then  x  is  just  the  fraction  of  the  total  post¬ 
collision  internal  energy  that  ends  up  in  the  first  molecule, 
x  is  distributed  with  a  probability  proportional  to  f(x), 
given  by 


f (x)  =  B  xc(l-x)d 


where 


B  =  [(<t>  ~  2) 


(<C>-2) 


]/[=c*d] 
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c 


Cx/2  -  1 


/ 


(46) 


and 

d  =  c2/  2  -  1  .  (47) 

Equation  (44)  can  be  sampled  via  the  acceptance-rejection  method 
to  determine  the  allocation  of  the  internal  energy  for  the  gen¬ 
eral  case.  For  the  special  case  of  both  molecules  being  diatomic, 
Eq.  (44)  becomes  singular,  with  the  limit  being  the  trivial  case 
that  all  distributions  of  internal  energy  are  equally  likely, 
i.e. , 

x  =  3  ,  (Ci  =  C2  =  2)  .  (48) 

For  the  special  case  of  just  one  of  the  molecules  being  diatomic 
(molecule  2,  for  instance),  then  the  distribution  for  x  can  be 
sampled  directly  via 

(2/Ci) 

x  =  8  1  ,  (C2  =  2)  ,  (49) 

with  the  obvious  reciprocal  relation  applying  for  C^  =  2.  The 
SSI  codes  recognize  these  special  cases  so  that  the  sampling 
can  be  expedited  when  possible,  while  retaining  the  full  gener¬ 
ality  of  Eqs.  (37)  and  (44)  when  required. 

4.3  Collisions  for  Molecules  with  Distinct  Weighting  Factors 

There  is  an  obvious  problem  when  considering  a  collision 
between  two  simulated  molecules  with  distinct  weighting  factors, 
since  they  represent  a  different  number  of  real  molecules.  If 
Wy  and  WL  represent  the  weighting  factors  for  the  two  molecules, 
with  W  being  greater  than  WT ,  then  the  collision  is  always 
counted  as  WL  "events".  This  is  accomplished  by  always  assign¬ 
ing  post-collision  velocity  and  energy  components  to  the  state 
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vector  of  the  molecule  with  the  smaller  weighting  factor,  but 
only  changing  the  components  of  the  molecule  with  the  greater 
weighting  factor  some  of  the  time.  The  probability  that  the 
molecule  with  the  greater  weighting  factor  will  have  its  com¬ 
ponents  changed  is  simply  WL/Wy.  Statistically,  this  means 
that  for  a  large  number  of  simulated  collisions,  each  such 
simulated  collision  will  average  out  to  WL  real  collisions  for 
each  species,  even  though  their  weighting  factors  differ.  It 
should  be  noted  that  this  does  violate  conservation  of  momentum 
and  energy  on  an  individual  collision  basis,  but  these  quanti¬ 
ties  are  conserved  in  the  aggregate  over  a  large  number  of 
collisions. 


4.4  Reactive  Collisions 

A  realistic  simulation  of  chemical  reactions  is  a  crucial 
element  of  many  problems.  If  a  bimolecular  reaction  of  the 
form 

A  +  B  C  +  D  (50) 

has  an  Arrhenius  rate  constant  of  the  form 


k  =  ArTn  exp(-Ea/RQT) 


(51) 


then  it  is  possible  to  define  a  reactive  cross  section  as  a 
function  of  translational  collision  energy  such  that  the  above 
rate  constant  is  implied  for  a  gas  in  translational  equilibrium. 
(In  Eq.  (51) ,  T  is  temperature,  E=  is  the  activation  energy, 
n  is  a  dimensionless  exponent  and  Ar  is  a  prefactor.)  The 
product  of  reactive  cross  section,  a*,  times  relative  velocity 

3 

can  be  expressed 


vfo* 


(1  +  6^)  VTT  Af 
2Rq  T  (n  +  3/2) 


(52) 
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where  6^  is  unity  for  like  reactants  and  zero  for  unlike  reac¬ 
tants.  The  probability  of  a  collision  resulting  in  a  reaction 
is  given  by  the  ratio  of  the  reactive  to  the  gas  kinetic  cross 
section  at  the  existing  relative  velocity  between  the  molecules 

In  the  SSI  Monte  Carlo  codes,  the  following  procedure  is 
used  to  determine  which,  if  any,  reaction  will  occur: 

•  All  possible  reactions  are  identified  for  the  pair 
of  molecules  which  are  selected  to  experience  a 
collision.  If  there  are  no  such  reactions,  the 
procedure  is  bypassed. 

•  The  probability  of  each  allowed  reaction  is  calcu¬ 
lated  by  computing  the  reactive  to  gas  kinetic 
cross  section  ratio. 

•  If  somehow  the  reactive  cross  section (s)  total 
to  a  greater  value  than  the  gas  kinetic  cross 
section  (which  will  rarely  be  the  case)  then 
the  probabilities  are  normalized  by  their  sum. 

•  The  collision  is  taken  to  be  reactive  with  a 
probability  equal  to  the  sum  of  the  individual 
reaction  probabilities. 

•  If  the  collision  is  reactive,  the  reaction  that 
occurs  is  selected  in  accordance  with  its 
probability. 

•  If  no  reaction  is  decided  upon,  then  the  colli¬ 
sion  is  either  elastic  or  inelastic  according 
to  the  user  specified  probability. 

Equation  (52)  becomes  singular  for  n  <  -3/2,  and  this 
procedure  cannot  handle  that  case.  The  major  limitation  of 
this  procedure  is  that  it  ignores  the  effect  of  internal  energy 


in  determining  whether  or  not  a  reaction  can  occur,  assuming 
that  the  entire  activation  energy  barrier  must  be  overcome 
through  translational  energy.  If  a  specified  fraction  of  the 
molecular  internal  energy  is  allowed  to  contribute  to  overcoming 
this  barrier,  then  the  restriction  on  n  is  relaxed  somewhat  (see 
Ref.  3),  but  this  feature  is  not  presently  implemented  in  the 
SSI  codes. 

4.5  Post-Reaction  State  Definition 

The  post-reaction  state  is  simply  determined  by  adding  the 
heat  of  reaction  to  the  total  collision  energy  (translational 
plus  internal)  and  then  computing  reactant  state  vectors  for 
the  products  based  on  this  total  energy  and  the  reactants 
center-of-mass  velocity,  as  for  inelastic  collisions.  Although 
this  is  clearly  an  oversimplification,  it  is  quite  consistent 
with  the  phenomenological  nature  of  the  model. 

The  position  state  vector  components  for  the  products  are 
randomly  selected  from  the  position  state  vector  components  of 
the  reactants,  which  are  not  the  same  (see  Section  6).  The 
major  additional  complication  of  chemical  reactions  is  that  of 
distinct  weighting  factors.  Since  the  reaction  is  V»L  "events" 
(see  Section  4.3),  it  only  destroys  the  reactant  with  the 
greater  weighting  factor  with  a  probability  of  WL/Wu*  If  Wp 
is  a  product  weighting  factor,  it  is  necessary  to  produce  WT/W 

h  p 

simulated  molecules  of  that  product.  In  general  this  is  not 
an  integer  quantity,  so  it  is  necessary  to  interpret  the  ratio 
statistically  so  that  the  expectation  value  is  proper.  That  is, 
sometimes  the  next  lower  integer  is  selected  and  sometimes  the 
next  higher  one,  with  a  probability  that  reflects  how  close  each 
integer  value  is  to  the  desired  fractional  quantity.  (See  the 
discussion  of  molecular  cloning  in  Section  5.2.)  Of  course,  the 
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weighting  factors  for  the  two  products  will  in  general  be  dif¬ 
ferent,  resulting  in  a  different  number  of  simulated  molecules 
being  produced  for  the  two  products.  Sometimes  this  process 
could  result  in  a  very  large  production  of  simulated  molecules 
for  a  product  with  a  small  weighting  factor.  In  order  to  pre¬ 
vent  overflow  of  code  dimensions  it  is  necessary  for  the  code 
to  sense  when  this  is  happening  and  automatically  increase  the 
product  weighting  factor  to  prevent  it. 

4.6  Dissociative  Reactions 

The  situation  for  dissociative  reactions  is  somewhat  com¬ 
plicated  by  the  presence  of  three  rather  than  two  products.  With 
a  little  manipulation,  however,  it  is  possible  to  use  the  previous 
relations  for  this  case  as  well.  Let  M^,  M2  and  represent  the 
three  product  molecules  from  a  dissociative  reaction,  with  a 
known  center-of-mass  velocity  and  total  energy.  The  procedure 
for  defining  the  post-reaction  state  is  to  first  define  an  arti¬ 
ficial  complex  comprised  of  the  (M2,M3)  pair.  (There  is  no 
implication  that  there  actually  is  any  such  collision  complex.) 

The  complex  is  assigned  a  number  of  internal  degrees  of  freedom 
equal  to  £  ,  given  by 

C 

cc  =  c2  +  C3  +  2(2  -  u>)  ,  (53) 

where  the  last  term  represents  the  contribution  of  the  relative 
translation  between  M2  and  M^. 

The  fact  that  this  term  is  not  simply  three,  as  it  would 
seem  it  should  be,  merits  some  explanation.  It  is  due  to  the 
fact  that  these  molecules  are  not  random  samples  from  the  gas 
but  rather  special  molecules  owing  to  their  being  created  in  a 
reaction.  This  point  can  perhaps  best  be  seen  by  considering 
microscopic  reversibility,  where  the  inverse  reaction  is  a 


30 


three  body  recombination.  For  this  reverse  process,  molecules 
participating  in  it  are  not  all  equally  probable,  since  those 
with  greater  relative  velocities  are  more  likely  to  collide. 
Hence,  the  term  does  take  on  the  value  three  for  the  special 
case  of  a)  equal  to  1/2,  which  is  precisely  the  case  of  collision 
frequency  being  independent  of  relative  velocity.  Examination 
of  Eqs.  (37)  and  (44)  demonstrates  that  translational  energy 
in  collisions  behaves  like  another  source  of  internal  energy 
with  2(2  -  w)  degrees  of  freedom. 

With  the  number  of  degrees  of  freedom  defined,  the  separa¬ 
tion  of  M^  from  the  (M2,M3)  complex  is  treated  as  an  inelastic 
collision.  The  resulting  velocity  for  the  complex  is  then 
treated  as  the  center-of-mass  velocity  for  M2  and  M^,  and  the 
internal  energy  assigned  to  the  complex  becomes  the  total  energy 
for  the  pair.  Using  these  values,  M2  and  M3  are  then  separated 
in  another  application  of  the  rules  for  inelastic  collisions. 

* 

5.  MOLECULAR  TRANSLATIONS 

5.1  Translation  in  an  Axisymmetric  Coordinate  System 

As  discussed  in  Section  1,  an  essential  element  of  the 
direct  simulation  Monte  Carlo  method  is  the  periodic  advancement 
of  simulated  molecules  along  their  trajectories.  Formally,  this 
is  accomplished  by  updating  the  position  and  velocity  elements 
of  the  state  vector.  For  a  cartesian  coordinate  system  this  is 
a  trivial  process,  but  the  relations  are  slightly  more  compli¬ 
cated  for  the  often  used  axisymmetric  coordinate  system.  Let 
vr0'  %o'  and  Vz0  rePresent  t*ie  initial  radial,  azimuthal  and 
axial  velocity  components  of  a  molecule,  with  r^  and  Zq  repre¬ 
senting  its  initial  radial  and  axial  position.  Additionally, 
let  $q  represent  the  initial  azimuthal  angle  for  the  molecule. 


This  is  included  here  just  for  demonstration  purposes;  it  will 
not  generally  be  known,  nor,  as  will  be  shown,  will  it  be 
needed . 


The  initial  position  and  velocity  of  the  molecule  can  then 
be  referenced  to  a  standard  cartesian  coordinate  system,  yielding 

?0  ■  [vr0  cos<*0>  •  v*0  sin(VJ* 

+  [vr0  sln(*0)  +  v0O  cos<y]j 
+  vz0  * 

*  vx0  1  +  vy0  i  *  vz0  k  1541 

and 

r0  "  [ro  c°s(<l>0)Ji  +  [r0  sin(<^0)Jj  +  zQ  k 

85  x0  1  +  y0  j  +  z0  *  *  (55) 


After  a  time  increment  Dt,  the  position  vector  of  the  molecule 
will  be 

*1  =  (x0  +  vx0  Dt)  *  +  &0  +  vy0  Dt)  ?  +  (20  +  vz0  Dt)  * 

-k  -k 

**  xi  i  +  yj  j  +  2i  k  ,  (56) 

and,  in  the  absence  of  perturbing  forces,  the  velocity  vector 
will  remain  unchanged  in  the  cartesian  coordinate  system.  It 
will  change  in  the  axisymmetric  coordinate  system  since  the 
basis  vectors  are  a  function  of  position  and  the  molecule  has 
moved.  The  new  radial  position  is 


(r0  +  vrO  Dt)  +  (v*0  Dt) 


and  the  new  axial  position  is 


*  *0  +  vz0  Dt 


The  radial  velocity  component  in  the  coordinate  system  appro¬ 
priate  to  the  new  molecular  position  can  be  determined  by  noting 


vrl  "  (v0>  '  (rl)/rl 


■  [ 
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vr0  (r0  +  vr0Dt)  +  %0(v»0Dt) 


and  it  can  similarly  be  shown  that 

V  -  v*0<Vrl> 


vzl  "  vz0  •  (61) 

Equations  (57)  -  (61)  give  the  updated  position  and  velocity 
elements  of  the  state  vector  for  a  translation  in  an  axisymmetric 
coordinate  system.  Note  that  these  relations  are  indeed  independ 
ent  of  <Pq.  A  similar  procedure  applies  for  molecular  transla¬ 
tion  in  other  coordinate  systems. 


5.2  Molecular  Cloning 

When  a  simulated  molecule  is  translated  from  one  cell  to 
another,  the  weighting  factor  for  that  species  will  generally 
be  different  in  the  new  cell.  Since  it  is  the  number  of  real 
molecules  rather  than  the  number  of  simulated  molecules  which 


mast  be  preserved  when  crossing  cell  boundaries  (statistically, 
at  least),  it  is  necessary  to  correct  for  the  distinct  weighting 
factors  (see  Section  3.4). 

If  the  weighting  factor  before  translation  is  Wq,  then  the 
simulated  molecule  represents  that  many  real  molecules.  If  the 
weighting  factor  in  the  new  cell  is  W^,  then  W^/W^  simulated 
molecules  would  be  required  to  represent  the  same  number  of 
real  molecules  in  the  new  cell.  If  this  ratio  were  a  whole 
integer,  then  this  could  be  accomplished  by  introducing  that 
many  "clones"  of  the  molecule  in  the  new  cell.  That  is,  that 
many  simulated  molecules  would  be  placed  in  the  new  cell,  all 
with  the  same  state  vector. 

When  the  number  WqA^  is  not  an  integer  (the  usual  case, 
of  course),  then  the  cloning  must  be  done  on  a  statistical 
basis.  So,  for  instance,  if  WQ/W1  were  equal  to  2.7,  then  30% 
of  the  time  two  clones  would  be  produced  and  70%  of  the  time 
three  would  be  produced.  Note  that  the  ratio  may  be  less  than 
unity,  and  the  molecule  may  not  be  introduced  into  the  new  cell 
at  all.  (In  which  case  the  molecule  is  removed  from  the 
simulation. ) 

Cloning  is  a  necessary  evil  inherent  in  a  system  with  spa¬ 
tially  dependent  weighting  factors.  It  enables  such  a  system 
to  maintain  the  statistically  correct  flux  of  mass  and  momentum 
across  cell  boundaries,  but  it  misrepresents  the  flux  of  random¬ 
ized  or  thermal  energy.  This  can  be  seen  by  an  extreme  case 
where  a  very  large  number  of  clones  is  produced  when  a  simulated 
molecule  crosses  a  cell  boundary.  The  resulting  molecules  in 
the  new  cell  have  the  correct  mass  and  momentum  flux,  but  since 
they  all  have  precisely  the  same  velocity  they  have  a  null 
relative  velocity,  and  therefore  a  zero  temperature.  If  the 
weighting  factors  are  not  too  different  between  adjacent  cells. 


then  the  errors  introduced  by  this  process  are  acceptably  small. 
However,  it  does  mean  that  one  cannot  arbitrarily  improve 
statistics  in  one  portion  of  the  solution  region  by  selectively 
reducing  the  weighting  factors  there.  This  was  a  difficulty 
which  was  encountered  in  the  early  stages  of  the  direct  simu¬ 
lation  Monte  Carlo  method  while  trying  to  improve  statistics 
along  the  axis  of  axisymmetric  simulations,  since  the  cell 
volumes  (and  therefore  the  sample  sizes)  tend  to  be  smallest 
on  the  axis. 

As  was  the  case  for  simulated  molecules  produced  via 
chemical  reactions,  it  is  possible  for  the  weighting  factors 
between  successive  cells  to  be  so  different  that  a  prohibitively 
large  number  of  simulated  molecules  would  be  required  to  produce 
the  same  number  of  real  molecules.  This  is  most  frequently  the 
case  when  a  new  species  is  being  introduced,  since  before  the 
species  gets  there  the  weighting  factor  is  initialized  to  a 
very  small  number.  The  codes  sense  when  a  disproportionate 
number  of  simulated  molecules  are  being  produced  for  a  given 
species  and  cell,  and  adjust  the  weighting  factor  automatically. 
As  the  weighting  factor  is  increased,  a  proportionate  fraction 
of  molecules  of  that  species  and  cell  are  removed  from  the 
simulation  in  order  to  keep  the  number  of  real  molecules  properly 
represented.  This  process  enables  the  weighting  factors  to  seek 
their  own  proper  level  without  a  priori  knowledge  of  the  solu¬ 
tion.  (Periodically,  the  cells  are  examined  to  determine  if  a 
certain  species  has  been  underrepresented  in  terms  of  number  of 
simulated  molecules.  If  it  is  found  to  be  the  case,  then  the 
weighting  factor  is  decreased,  allowing  weighting  factors  to 
float  downwards  as  well  as  upwards.  It  is  the  danger  of  weight¬ 
ing  factors  b  .mg  too  small,  causing  an  overflow  of  code  dimen¬ 
sions,  which  is  most  critical,  however.) 

i 
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COLLISION  SAMPLING  IN  A  MULTI -COMPONENT  VHS  GAS 


6.1  General  Considerations  and  Approach 

The  two  general  considerations  in  the  sampling  of  colli¬ 
sions  are,  as  usual,  accuracy  and  efficiency  of  the  simulation. 
As  far  as  accuracy  is  concerned,  it  is  crucial  that  the  method 
in  which  molecules  are  selected  for  collisions  be  proper.  It 
is  imperative  that  the  correct  collision  frequencies  be  simu¬ 
lated  between  various  species,  and,  in  fact,  between  the  dif¬ 
ferent  portions  of  the  velocity  phase  space  for  the  various 
species.  Furthermore,  this  frequency  of  simulated  collisions 
must  remain  correct  without  any  requirements  put  on  the  velocity 
distribution  function;  it  certainly  must  not  be  assumed  that 
there  is  a  Maxwellian  velocity  distribution. 

As  far  as  efficiency  is  concerned,  it  is  highly  desirable 
to  use  a  method  of  collision  sampling  involving  a  computational 
effort  which  is  proportional  to  the  number  of  simulated  mole¬ 
cules,  N,  in  a  cell.  Methods  which  are  proportional  to  a  power 
of  N  greater  than  unity  can  become  prohibitively  time  consuming 
as  the  number  of  molecules  is  increased  —  a  limit  which  should 
be  made  as  accessible  as  possible  for  obvious  physical  reasons. 

6.2  Collision  Sampling  for  a  Single  Component  Gas 

The  simplified  situation  of  a  simulation  involving  only 
one  species  is  considered  here.  This  problem  is  significant  in 
part  due  to  all  the  attention  it  has  received  and,  as  will  be 
seen,  it  serves  as  an  important  reference  case.  When  there  is 
just  one  species,  then  there  is  just  one  gas  Kinetic  cross 
section  (though  it  is  still,  of  course,  a  function  of  collision 
energy),  just  one  molecular  weight  and  just  one  weighting  factor 
for  each  cell.  In  short,  just  one  of  everything  that  has  a 


molecular  subscript.  Hence,  in  this  subsection  all  such  quanti¬ 
ties  will  be  presented  without  subscripts.  The  most  important 
simplification  of  having  a  single  species  is  that  there  is  just 
one  collision  class,  i.e.,  only  self-collisions  of  the  given 
species  with  itself  are  possible. 

6.2.1  Collision  Pair  Selection 

As  discussed  in  Section  1,  in  the  direct  simulation  Monte 
Carlo  method  collisions  are  sampled  on  a  cell  by  cell  basis 
until  the  number  of  collisions  simulated  is  appropriate  to  the 
overall  solution  time  step,  Atm»  The  only  spatial  requirement 
placed  on  potential  collision  partners  is  that  they  be  within 
the  same  cell.  In  particular,  it  is  not  required  that  they  be 
within  a  molecular  diameter  of  each  other.  (Note  that  if  all 
pairs  of  molecules  were  inspected  to  find  those  that  were  suf¬ 
ficiently  close  to  each  other,  this  would  involve  a  computa¬ 
tional  effort  in  proportion  to  the  square  of  the  number  of 
molecules  in  the  cell.)  The  rationale  for  this  is  that  the 
cells  should  be  taken  small  enough  such  that  macroscopic  prop¬ 
erties  can  be  assumed  constant  across  the  cell.  When  this  is 
the  case,  then  a  molecule  within  the  cell  can  be  considered 
typical  of  a  molecule  which  might  exist  anywhere  within  the 
cell,  and  molecular  location  can  be  ignored  when  selecting 
potential  collision  pairs. 

Spatial  consideration  aside,  the  probability  of  any  two 
molecules  experiencing  a  collision  is  proportional  to  acr,  the 
product  of  their  mutual  cross  section  times  their  relative 
velocity.  This  probability  is  correctly  simulated  via  an 
application  of  the  acceptance-rejection  method  if  pairs  of 
molecules  are  selected  at  random,  and  then  kept  or  rejected  as 
collision  partners  with  a  probability  proportional  ocr.  This 


is  accomplished  by  keeping  a  maximum  value  for  ocr  for  each 
cell  (which  is  updated  if  a  larger  value  is  encountered)  and 
computing  the  ratio  r  defined  by 


A  random  variable,  3,  is  then  generated,  and  the  pair  of  mole¬ 
cules  is  accepted  as  collision  partners  if  r  is  greater  than  3. 
This  produces  the  proper  relative  collision  probability  without 
regard  to  the  existing  velocity  distribution  function. 

6.2.2  Collision  Time  Counter  for  a  Single  Component  Gas 

The  volumetric  collision  frequency  for  a  single  component 
gas,  v  (collisions  per  unit  volume  per  unit  time),  is  given  by 

1  2  _ 

v  —  2"  n  ocr  •  (63) 

where,  as  in  Section  2,  n  represents  the  number  density  of  the 
species  and  ocr  is  the  average  product  of  collision  cross  sec¬ 
tion  and  relative  velocity.  At  first  inspection,  it  would  seem 
from  Eq.  (63)  that  a  correct  simulation  of  collision  frequency 
would  require  evaluation  of  ocr ,  which  would  mean  that  all  pairs 
of  molecules  in  a  cell  would  have  to  be  considered.  Such  a 

2 

procedure  involves  a  computational  effort  proportional  to  N 
and  is  to  be  avoided,  if  possible,  in  preference  to  a  method 
which  is  simply  proportional  to  N. 

The  alternative  approach,  introduced  by  Bird,  is  the  time 
counter  approach.  For  each  collision  a  time  counter,  tc,  is 
incremented  by  an  amount  that  depends  on  the  relative  velocity 
of  the  collision.  Collision  sampling  continues  in  a  cell  until 
its  time  counter  has  been  advanced  beyond  the  overall  flow 
simulation  time,  at  which  time  the  code  proceeds  to  the  next 


cell  (which  has  its  own  time  counter) .  The  time  counter  incre¬ 
ment,  At  is  given  by 

v 


At  = 


2 

Vn  crc 


(64) 


where  V  is  the  cell  volume  and  n  is  the  species  number  density 
given  by 

n  =  NW/V  ,  (65) 

with  W  being  the  weighting  factor  for  the  species.  (Equation 
(65)  is  just  a  special  case  of  Eq.  (26).)  It  should  be  stressed 
that  Eq.  (64)  applies  for  each  real  collision.  As  is  discussed 
in  Sections  3.4  and  4.3,  each  simulated  collision  corresponds  to 
W  real  collisions,  so  when  a  simulated  collision  occurs  the 
actual  applied  increment  to  t  is  W  times  the  value  given  by 
Eq.  (64). 

It  is  not  obvious  that  Eq.  (64)  will  lead  to  a  proper 
simulation  of  the  overall  collision  frequency,  so  a  demonstra¬ 
tion  will  be  presented.  Let  f-^tc^.)  be  the  normalized  distribu¬ 
tion  function  for  relative  velocity  appropriate  to  a  given  cell 
at  a  given  time.  (Note  that  most  problems  solved  by  a  Monte 
Carlo  technique  involve  repeated  runs,  where  the  total  number 
of  collision  pairs  can  be  made  arbitrarily  large.  The  intro¬ 
duction  of  a  distribution  function,  and  therefore  the  demonstra¬ 
tion  of  the  correctness  of  Eq.  (64)  is  strictly  valid  only  in 
the  limit  of  infinitely  many  runs.  Since  the  essence  of  the 
method  is  a  correct  statistical  representation,  it  could  not 
be  otherwise.)  By  definition, 


J  fl(cr)dcr 


1 


(66) 
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irrespective  of  the  particular  form  for  f^  Let  f2(cr)  be  the 
normalized  distribution  function  of  relative  velocity  occurring 
in  collisions.  Since  collisions  occur  with  a  probability  pro¬ 
portional  to  ocr,  f 2  can  be  constructed  from  f^  via 


f2(cr)  =  fl(cr)  ^ 

oc 


where  acr  can  now  be  formally  defined  via 


=r  "  /  fl(cr>ocrdcr 


(Note  that  in  Eq.  (68),  as  in  the  rest  of  this  demonstration, 
the  functional  form  of  the  dependence  of  cross  section  on  rela¬ 
tive  velocity  need  never  be  specified.  The  time  counter  repre¬ 
sented  by  Eq.  (64)  is  therefore  not  restricted  to  any  particular 
model  for  the  cross  section.) 

The  average  increment  of  the  time  counter  that  is  applied 
over  many  simulated  collisions  is  therefore  given  by 


"  /  f2(cr)Atc(cr)dcr  * 


If  Eqs.  (64)  and  (67)  are  substituted  into  Eq.  (69),  the  result 


n^Voc 


where  the  normalization  condition  (Eq.  (66))  has  been  utilized. 
The  implication  of  Eq.  (70)  is  that,  on  average,  the  frequency 
of  collisions  in  the  cell  is  Vn^ocZ/2.  Since  this  is  simply 


the  product  of  cell  volume  times  the  proper  volumetric  colli¬ 
sion  frequency  (Eq.  (63)),  the  validity  of  the  time  counter 
given  in  Eq.  (64)  has  been  demonstrated. 


6.3  Collision  Class  Sampling  in  Gas  Mixtures 

The  above  procedure  for  a  single  species  gas  can  be  extended 
to  a  multi-component  mixture  via  consideration  of  distinct  col¬ 
lision  classes.  In  this  approach,  collision  classes  are  defined 
by  the  colliding  pair  identities.  Hence,  if  there  are  p  species 
in  the  simulation  then  there  are  p(p+l)/2  collision  classes, 
which  can  be  identified  by  the  subscripts  of  the  corresponding 
molecular  pair.  (The  number  of  classes  is  not  p  since  the 
order  of  molecule  specification  is  not  taken  to  matter  in  deter¬ 
mining  a  collision  class.  Hence,  the  class  identified  by  the 
subscripts  i,j  is  not  distinct  from  the  class  identified  by  the 
subscripts,  j,i.) 

In  collision  class  sampling,  which  is  the  method  used  by 

Bird,"*"  each  collision  class  is  sampled  separately,  and  the 

collision  sampling  in  a  cell  is  not  complete  until  all  classes 

have  been  considered.  Each  collision  class  has  its  own  stored 

value  of  (a- .c  )  _ ,  and  its  own  separate  time  counter,  t  . .. 

rj  r  max  ci} 

By  a  comparable  analysis  to  that  presented  above,  it  can  be 
shown  that  the  appropriate  time  counter  increment  in  this  case 
is 


a  ♦  «„) 


ninjVoi jcr 


where,  as  in  Section  4,  is  the  Kronecker  delta  which  is 
unity  for  i  =  j  and  zero  otherwise.  As  in  the  previous  section, 
the  above  increment  applies  for  each  real  collision.  A  simu¬ 
lated  collision  corresponds  to  WT  real  collisions,  where  W.  is 

Li  J-* 
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the  lesser  of  and  (see  Section  4.3),  so  when  a  simulated 

collision  occurs,  the  applied  increment  to  t  .  t  is  WT  times  the 

ci  ^  L 

result  of  Eq.  (71) . 

6.4  Global  Collision  Sampling  in  a  Gas  Mixture 

Although  the  procedure  described  above  is  quite  reasonable 
for,  say,  a  two  component  mixture,  it  becomes  quite  complicated 
as  the  number  of  species  increases.  For  10  species,  for  instance, 
the  program  must  loop  over  55  distinct  collision  classes  for  each 
cell,  and  storage  must  be  allocated  for  110  quantities  in  each 
cell.  As  the  number  of  species  increases,  the  storage  require¬ 
ment  for  the  collision  sampling  constants  quickly  becomes 
greater  than  the  storage  required  for  the  molecular  state  vectors! 
The  obvious  simplification  is  to  search  for  a  technique  where 
collisions  are  simulated  simultaneously  for  all  collision  classes, 
with  each  class  having  its  proper  relative  probability  of  being 
selected.  The  overall  collision  sampling  then  continues  until  a 
single  time  counter  indicates  that  sufficient  collisions  have 
been  sampled  in  the  current  time  step  and  cell. 

6.4.1  Global  Collision  Time  Counter 

If  molecular  pairs  are  selected  for  collisions  such  that 
the  various  collision  classes  automatically  appear  with  the 
proper  relative  frequency  (see  below),  then  it  is  not  necessary 
to  consider  separate  time  counters  for  all  the  various  collision 
classes.  One  approach  that  could  then  be  applied  is  to  keep  a 
collision  time  counter  for  just  one  collision  class,  and  incre¬ 
ment  it  when  collisions  of  that  class  occur.  If  the  various 
collision  classes  are  being  selected  according  to  their  correct 
relative  frequency,  then  simulating  the  proper  frequency  for 
one  collision  class  will  ensure,  in  the  long  run,  that  all 


collision  classes  are  occurring  with  the  correct  frequency.  A 
disadvantage  with  this  approach  is  the  necessity  of  making  an 
arbitrary  choice  for  the  collision  class  which  is  to  have  a 
time  counter.  Furthermore,  there  may  be  no  good  choice  in  a 
reacting  flow  where  the  dominant  species  may  vary  strongly 
from  place  to  place.  (Clearly,  one  would  not  want  to  select  a 
class  of  collision  that  does  not  occur  in  a  given  cell,  since 
the  result  would  be  a  never  ending  sampling  of  collisions  of 
other  classes.) 

The  preferred  approach  is  to  define  a  global  collision 
time  counter,  tg,  which  is  a  weighted  average  of  the  time 
counters  of  all  collision  classes;  i.e. 


P  i 


E  E° 

W . 

p  1 


.  .  t  .  . 
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E  E°ij 

i=l  j=l 

where  the  are  nonnegative  coefficients  which  can  be  selected 
at  will.  Note  that  in  this  formulation  every  collision  will 
result  in  some  increment  of  the  global  time  counter  (unless 
Dij  *  ®  ^or  t^iat  class),  so  the  collision  sampling  frequency 
is  not  dependent  on  any  one  collision  class. 

It  remains,  of  course,  to  specify  the  ^ .  A  very  conven¬ 
ient  choice  is  given  by 


ninj 


Firstly,  Eq.  (73)  is  convenient  because  it  tends  to  make  the 
collision  classes  with  the  higher  collision  frequencies  count 
more,  resulting  in  good  statistics  for  t  irrespective  of  cell 
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location.  (Note  that  ^  is  cell  dependent  since  the  species 
number  densities  are  cell  dependent.)  Secondly,  Eq.  (73) 


results  in  a  particularly  convenient  form  for  t^.  The  normal¬ 


izing  factor  in  the  denominator  of  Eq.  (72)  can  be  summed  to 
give 
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T  fccij 


(74) 


Hence,  a  collision  of  class  i j ,  which  would  produce  an  incre¬ 


ment  of  At  . .  to  its  own  time  counter  produces  an  increment  At 
ci3  r  g 

to  t  given  by 


At 


At 


'ci  3 


(75) 


where,  again,  n  is  the  total  number  density  of  all  species  in 
the  cell.  If  eq.  (71)  is  substituted  into  Eq.  (75),  the  result 
is 


At. 


(76) 


Vn  °ijcr 


Equation  (76)  is  extremely  significant  since  it  recaptures  the 
precise  form  of  the  time  counter  increment  for  a  single  species 
(Eq.  (64)),  but  indicates  that  it  is  completely  valid  for  a 
multi-component  mixture  so  long  as  the  various  collision  classes 
are  sampled  with  the  proper  relative  frequency. 


6.4.2  Collision  Pair  Selection  in  Multi-Component  Mixtures 

When  considering  selection  of  collision  pairs,  it  is  cru¬ 
cial  to  remember  the  distinction  between  real  and  simulated 
molecules  discussed  in  Section  3.4.  Given  two  simulated 


molecules  selected  at  random  from  within  the  cell,  the  probabil¬ 
ity  of  their  having  a  real  collision  is  proportional  to  WiWj°ijcr* 
However,  real  collisions  cannot  happen  individually;  they  come 
W_  at  a  time,  where  WT  is  the  lesser  of  W.  and  W..  Hence,  when 

Lt  Li  X  j 

a  collision  is  decided  upon  in  the  program,  WL  of  them  will 
occur.  To  compensate  for  this,  potential  collision  pairs  should 
be  accepted  for  collision  according  to  the  size  of  Q  given  by 


Q 


Vijcr 


(77) 


The  relative  frequency  of  real  ij  collisions  will  then  be  pro¬ 
portional  to  the  product  QWL  (the  relative  probability  of  a  pair 
being  accepted  for  collision  times  the  number  of  real  collisions 
occurring  when  the  pair  is  accepted),  which  is  the  desired  rela¬ 
tion.  Selection  of  collision  pairs  with  the  correct  relative 
frequency  then  assures  that  incrementing  the  global  time  counter 
as  discussed  above  will  give  a  statistically  correct  sampling 
of  all  collision  classes  simultaneously. 


6.4.3  Summary  of  Collision  Sampling  in  Multi-Component 
Mixtures 

The  results  of  this  section  can  be  summarized  via  the  fol¬ 
lowing  procedure  for  the  sampling  of  collisions: 

•  Each  cell  has  a  (current)  maximum  value  of  Q, 

Qmax'  that  has  been  encountered  so  far  in  the 
collision  sampling  process.  Whenever  a  larger 
value  is  encountered,  Qmax  is  set  equal  to  that 
larger  value. 

•  Each  cell  has  a  current  value  of  the  global 
time  counter,  tg. 

•  Pairs  are  selected  at  random  from  all  molecules 
within  the  cell. 


•  For  each  pair,  Q  vas  defined  by  Eq.  (77))  is 
computed . 


•  The  ratio  of  Q  to  Qjnax  is  computed,  and  a  random 
variable  is  generated.  The  pair  is  accepted  for 
collision  if  the  random  variable  is  less  than 
the  ratio.  (If  the  pair  is  not  accepted,  then 
another  random  pair  is  selected.  The  process 
continues  until  a  pair  is  accepted.) 

•  For  an  accepted  pair,  the  collision  mechanics 
are  computed  as  described  in  Section  4.  This 
always  corresponds  to  Wl  collisions. 

•  The  global  time  counter  is  incremented  by 
WLAtg,  where  Atg  is  given  in  Eq.  (76). 

•  The  process  continues  until  the  global  time 
counter  goes  beyond  the  overall  flow  time. 

At  that  point,  the  collision  sampling  is  com¬ 
menced  in  the  next  cell. 

•  When  all  cells  have  had  collisions  simulated, 
then  the  code  proceeds  to  the  translation  portion. 
(See  Sections  1  and  5.) 


7.  INITIAL  AND  BOUNDARY  CONDITIONS 
7.1  General  Considerations 

The  initial  and  boundary  conditions  necessary  to  simulate 
a  problem  frequently  do  not  receive  their  fair  share  of  consid¬ 
eration.  It  is  these  conditions  which  usually  distinguish  one 
solution  from  another,  and  their  correct  and  efficient  specifi¬ 
cation  should  be  a  central  concern.  Nonetheless,  there  is 
clearly  room  for  advancements,  particularly  in  the  specification 
of  boundary  conditions.  Many  gas  dynamic  solutions  involve 
boundary  conditions  specified  at  infinity,  which  are  currently 
simulated  by  placing  boundaries  very  far  from  the  main  flow 
region.  It  would  result  in  a  substantial  computational  simpli¬ 
fication  if  boundary  conditions  applicable  closer  in  to  the  flow 
region  of  interest  could  be  generated  for  free  gas  boundaries. 
Wall  boundary  conditions  also  frequently  involve  a  fair  degree 


of  approximation ,  usually  taking  the  form  of  accommodation 
coefficients  (though  in  this  case  the  problem  is  as  much  a  lack 
of  basic  physical  understanding  of  the  gas-surface  interaction 
process  as  it  is  a  lack  of  a  good  numerical  simulation) . 

7.2  Initial  Conditions 

Since  the  direct  simulation  Monte  Carlo  method  is  inherently 
an  unsteady  technique  an  initial  state  must  be  specified  in 
order  to  advance  the  solution.  (For  situations  where  a  steady 
state  result  is  desired,  it  is  obtained  as  the  long  time  solu¬ 
tion  to  an  unsteady  problem.  In  this  case  the  initial  condi¬ 
tions  have  no  effect  on  the  eventual  solution,  but  they  may  well 
have  an  impact  on  the  speed  with  which  that  state  is  achieved.) 

It  will  be  assumed  here  that  the  initial  conditions  correspond 
to  a  uniform  flow  with  the  translational  and  internal  modes 
being  in  equilibrium.  The  specification  of  the  initial  condi¬ 
tions  therefore  involves  determining  the  state  vector  elements 
consistent  with  this  condition  for  the  desired  number  of 
molecules. 

7.2.1  Number  of  Simulated  Molecules  and  Weighting  Factors 

The  desired  number  of  simulated  molecules  of  each  species 
in  each  cell  (referred  to  here  as  M  )  will  usually  be  an  input 
quantity.  (Typically,  simulations  aim  for  a  total  number  of 
molecules  per  cell  in  the  neighborhood  of  twenty.)  Given  the 
initial  number  density  to  be  simulated  for  a  species,  n^, 

(which  will  have  been  converted  in  the  code  to  internal  dimen¬ 
sions  -  see  Section  3)  the  weighting  factor  for  a  species  in 
a  cell  can  be  dervied  from  an  application  of  Eq.  (26)  to  give 


where  V  is  the  cell  volume.  If  a  species  is  not  initially  pre¬ 
sent  in  a  cell,  then  the  weighting  factor  is  set  to  a  very  small 
but  positive  number,  typically  that  which  would  correspond  to 
an  initial  mole  fraction  of  one  part  per  million  or  so.  A 
weighting  factor  of  zero  would  cause  an  attempt  to  divide  by 
zero  when  molecules  of  that  species  move  into  the  cell,  but  it 
is  generally  best  to  start  the  weighting  factors  off  small. 

As  the  solution  proceeds,  the  weighting  factors  are  automatically 
adjusted,  but  the  adjustment  upward  is  more  direct  and  immediate 
than  the  adjustment  downward  (see  Sections  3.4  and  5.2). 

7.2.2  Initial  Positions 

The  initial  molecules  assigned  to  a  cell  should  have  an 
equal  probability  of  being  placed  in  any  volume  element  of  the 
cell.  The  rules  for  accomplishing  this  will  change  with  the 
coordinate  system  being  used,  but  they  are  readily  derivable 
from  the  general  principle.  As  an  example,  consider  an  axisym- 
metric  simulation  where  a  cell  is  bounded  by  the  radial  posi¬ 
tions  r.^  and  r2  (r^  <  r2),  and  the  axial  coordinates  z-^  and  z2 
(z^  <  z^).  Since  the  basic  volume  element  in  this  coordinate 
system  is  2irrdrdz  =  iTd(r2)dz,  the  volume  elements  will  be 
sampled  with  equal  probability  if  r  and  z  are  sampled  randomly. 
Hence,  a  molecule  can  be  assigned  axial  and  radial  positions  via 

r  =  >/r12  +  0(r22  -  r^)  (79) 

and 

z  =  z1  +  g(z2  -  z1)  .  (80) 


(Recall  that  every  time  the  symbol  8  occurs  it  implies  a  separ¬ 
ate  random  variable.  In  particular,  one  should  certainly  not 


use  the  same  random  variable  to  determine  radial  and  axial 
coordinates.  The  second  application  would  hardly  qualify  as 
“random" . ) 


7.2.3  Initial  Velocity  Components 

The  thermal  velocity  components  for  a  molecule  in  transla¬ 
tional  equilibrium  (neglecting,  for  the  moment,  any  mean  flow 
contribution)  should  be  selected  from  the  normalized  Maxwellian 
velocity  distribution,  fg(v),  given  by 

fQ  *  — —  exp  [-(av)2]  ,  (81) 


where 


m  is  the  species  molecular  weight,  RQ  is  the  universal  gas  con¬ 
stant  and  T^  is  the  temperature.  Equation  (81)  applies  for  each 
of  the  molecular  velocity  components,  and  must  be  sampled  three 
times  for  each  molecule  that  comprises  the  initial  state  of  the 
simulation.  A  method  for  directly  sampling  from  this  distribu¬ 
tion,  as  presented  in  Ref.  1,  is 

A.  =  2ttB  ,  (83) 


f-log(B)/a 


v  =  A2sin(A^)  .  (85) 

After  the  thermal  velocity  components  are  determined  for  each 
molecule,  then  any  mean  flow  velocity  is  simply  added  on.  The 
velocities  are  then  transformed  to  internal  units  (see  Section 


7.2.4  Initial  Internal  Energies 

The  only  remaining  element  of  the  state  vector  to  be  speci' 
fied  is  the  internal  energy.  Internal  energies  for  a  gas  in 
equilibrium  are  distributed  according  to  the  normalized  distri¬ 
bution  function  fj  given  by 

,?/2-l 

^i  =  F TUT)  exP  r  (86 

where  £  represents  the  number  of  internal  degrees  of  freedom 
for  the  species  in  question,  T  is  the  gamma  function  and  £  is  a 
dimensionless  internal  energy,  i.e. 

K  =  VR0T«  '  (87 

where  Ej  is  the  internal  energy.  In  general,  sampling  of  Eq. 
(86)  must  be  done  via  the  acceptance-rejection  method.  In  the 
present  SSI  codes  r,  is  restricted  to  being  greater  than  or 
equal  to  two  (or  the  trivial  case  of  £  equal  to  zero,  which 
just  gives  Ej  identically  equal  to  zero).  If  £  is  precisely 
equal  to  two,  then  a  direct  sampling  of  the  internal  energy 
is  possible  via 

£  =  -log  ( B)  .  (C  =  2)  (88 

In  the  general  case  of  c  >2,  it  proves  convenient  to  first 
introduce  the  transformation  s  =  vr.  s  is  then  distributed  in 
proportion  to  the  distribution  g(s)  given  by 

g  (s)  =  2sU“1)  exp (— s2 )  .  (89 

Since  g(s)  is  to  be  sampled  via  the  acceptance-rejection 
method,  it  is  first  necessary  to  determine  its  maximum  value, 
g„„„  Standard  calculus  serves  to  show  that  g„„„  occurs  at 
s  =  s*,  where 


The  sampling  of  Eq.  (89)  proceeds  as  follows: 


•  A  value  of  s  is  selected  randomly  from  the  inter¬ 
val  Sjnin  to  s^,  where 


s 

max 


s*  +  5 


(92 


and 


s  =  greater  of  {0,s*  -  5}  .  (93 

mm 

(Note  that  in  this  interval  g(s)  goes  from  its 
maximum  value  to  a  value  on  the  order  of  10“^° 
its  maximum  value.  The  transformation  from  Eg. 

(86)  to  Eq.  (89)  was  made  mainly  to  achieve  a 
probability  function  which  dies  off  extremely 
rapidly  away  from  its  maximum  value,  so  that 
very  little  error  is  associated  with  considering 
a  finite  subinterval  for  the  sampled  variable.) 

•  g(s)/gmax  *s  calculated  via  Eq.  (91). 

•  A  random  variable  is  generated,  and  the  value  of 
s  is  kept  if  the  random  variable  is  less  than 
g(s)/gmax*  Otherwise,  the  procedure  is  repeated 
until  a  value  of  s  is  accepted. 

•  When  a  value  of  s  is  selected,  then  the  internal 
energy  is  given  by  sZrqTjo. 


As  for  all  the  initial  conditions,  the  codes  will  automatically 
then  express  the  values  in  internal  dimensions  (see  Section 


7.3  Wall  Boundary  Conditions 

Boundary  conditions  in  the  direct  simulation  Monte  Carlo 
technique  are  applied  for  both  wall  and  free  gas  boundaries. 

In  a  wall  boundary  condition,  whenever  a  molecule  is  simulated 
to  strike  the  wall  some  action  must  be  taken.  Unless  the  wall 
is  a  condensing  boundary  the  molecule  will  be  reflected,  and 
one  of  several  boundary  conditions  can  be  selected.  The  easiest 
condition  is  that  of  a  specularly  reflecting  wall,  where  the 
molecule  simply  is  replaced  by  its  mirror  image  to  keep  it 
within  the  solution  region.  Another  frequently  used  condition 
is  that  of  a  perfectly  accommodating  wall.  In  this  condition, 
a  molecule  is  reemitted  from  the  wall  after  being  selected  from 
a  distribution  characteristic  of  the  wall  temperature  and  veloc¬ 
ity.  To  date,  wall  boundary  conditions  have  not  been  required 
in  the  SSI  codes,  and  they  will  not  be  discussed  further  here. 
They  are  discussed  at  some  length  in  Ref.  1,  including  the 
intermediate  case  of  partial  accommodation. 

7.4  Free  Gas  Boundary  Conditions 

In  many  cases,  the  boundary  condition  is  meant  to  simulate 
a  region  of  uniform  equilibrium  flow.  Molecules  leave  the  solu¬ 
tion  region  in  the  normal  course  of  their  trajectories,  and  they 
simply  disappear  from  the  simulation.  Molecules  are  introduced 
from  outside  the  boundary  as  selected  from  distributions  appro¬ 
priate  to  incoming  molecules  in  the  undisturbed  flow.  It  should 
be  stressed  that  this  is  not  the  same  as  simply  sampling  a 
Maxwellian  velocity  distribution,  since  it  is  the  molecular 
flux  across  the  boundary  which  must  be  correctly  simulated. 
Hence,  molecules  with  a  large  component  of  velocity  inward  from 
the  boundary  are  more  likely  to  be  selected  than  they  would  be 
in  choosing  molecules  appropriate  to  a  static  distribution  as 
was  done  for  the  initial  conditions  described  above. 


52 


7.4.1  Incoming  Number  Flux 

The  first  requirement  is  to  determine  the  number  of  mole¬ 
cules  which  should  be  introduced  across  the  boundary  during  the 

solution  time  step  At  ,  The  incoming  number  flux,  q,  (molecules 

.  111  1 
per  unit  area  per  unit  time)  can  be  expressed 

q  *  n^  |  exp  (-w^)/VtT  w[l  +  erf(w)]  j  /2a  ,  (94) 

where 

w  =  auwcos(0)  ,  (95) 

and  0  represents  the  angle  between  the  inward  surface  normal  and 
the  mean  flow,  which  has  a  magnitude  u^  (i.e.,  u^cosf©)  is  the 
inward  component  of  the  mean  flow  velocity) ,  n^  is  the  far  field 
number  density  of  the  species  of  interest  and  a  is  given  in 
Eq.  (82). 

Equation  (94)  must  be  applied  for  each  cell  on  the  boundary 
and  for  each  species  that  exists  in  the  ambient.  The  number  of 
simulated  molecules  to  be  introduced  into  the  cell,  Nb,  is 
given  by 

Nb  =  <jWw  -  (96> 

where  Ac  is  the  area  of  the  cell  along  the  boundary  and  W  is 
the  weighting  factor  for  the  species  and  cell  in  question. 

7.4.2  Incoming  Molecular  Velocity  Components 

A  coordinate  system  should  be  set  up  locally  at  the  bound¬ 
ary  such  that  one  direction  is  in  the  direction  of  the  inward 
normal  and  the  other  two  directions  are  perpendicular  to  it. 
Velocity  components  are  first  determined  in  terms  of  this  local 
coordinate  system,  and  then  transformed,  if  necessary,  to  the 


main  code  coordinate  system.  In  the  local  coordinate  system, 
the  velocity  components  parallel  to  the  surface  are  determined 
as  discussed  above  for  initial  state  molecules,  but  the  inward 
component  of  velocity  must  be  selected  in  proportion  to  the 
distribution  h(v)  given  by 

h(v)  =  (v  +  w)  [exp  -(av)2]  ,  (97) 

which  must  be  sampled  via  the  acceptance-rejection  method. 

(Only  positive  v  is  allowed,  of  course,  from  the  definition 
of  the  coordinate  system.) 

7.4.3  Incoming  Molecular  Position 

An  initial  entry  position  should  be  selected  for  the  mole¬ 
cule  such  that  the  flux  is  randomly  distributed.  Assuming  there 

is  no  variation  of  9  across  A  this  means  that  each  area  element 

c 

of  the  exposed  cell  area  should  have  an  equal  probability  of 
being  selected.  Once  the  initial  entry  position  is  selected, 
then  the  molecule  should  be  translated  a  random  fraction  of  At 

m 

along  its  trajectory  to  determine  its  actual  location.  All  of 
the  considerations  discussed  in  Section  5  with  regard  to  molecu¬ 
lar  translations  apply  to  this  translation  and,  in  particular, 
it  must  be  possible  to  dynamically  adjust  the  weighting  factor 
as  required. 


7.4.4  Incoming  Molecular  Internal  Energies 

The  internal  energies  are  selected  as  described  in 
Section  7.2.4. 
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STATISTICAL  SAMPLING  OF  OUTPUT 


8.1  General  Considerations 

It  is  safe  to  say  that  the  molecular  state  vectors  as  they 
exist  in  the  computer  do  not  comprise  the  desired  output  of  the 
procedure.  With  rare  exceptions,  it  is  usually  macroscopic 
quantities  such  as  temperature,  density,  mean  flow  velocity, 
etc.,  which  are  of  interest  —  not  the  microscopic  quantities 
represented  by  the  state  vector  components  of  an  individual 
simulated  molecule.  The  generation  of  the  desired  output 
requires  that  the  macroscopic  quantities  of  interest  be  repre¬ 
sented  in  terms  of  statistical  sums  of  the  available  microscopic 
quantities;  and  it  is  the  main  purpose  of  this  section  to  pre¬ 
sent  these  correspondences.  All  sums  are  kept  in  terms  of  "real 
molecules  and  events,  i.e.,  the  current  weighting  factors  are 
included  in  the  sums.  This  is  essential  since  the  weighting 
factor  determines  the  statistical  importance  of  a  given  molecule 
Since  the  weighting  factors  are  dynamically  and  unpredictably 
adjusted  as  the  solution  progresses,  it  would  not  be  possible 
to  go  back  and  add  in  the  effect  of  weighting  factors  a 
posteriori. 

In  general,  it  must  be  decided  ahead  of  time  exactly  what 
output  is  desired  from  the  code,  and  therefore  what  statistical 
sums  should  be  kept  to  generate  it.  There  is  a  vast  amount  of 
potential  information  in  the  simulation,  and  it  is  not  reason¬ 
able  to  store  all  possibly  interesting  quantities  in  all  runs. 

On  the  other  hand,  it  is  wasteful  to  completely  rerun  a  case 
just  because  the  user  decids  there  was  an  additional  quantity 
he  was  interested  in.  The  selection  of  output  for  a  given  run, 
therefore,  unavoidably  requires  user  judgement.  Once  the  user 
has  decided  upon  the  required  output,  the  determination  of  which 


statistical  sums  should  be  kept  can  be  done  automatically  by 
the  code.  Care  is  taken  to  make  sure  that  a  statistical  sum 
is  not  duplicated  internally  if  it  is  required  by  more  than  one 
requested  output  quantity. 

Some  initial  words  of  caution  are  required.  By  its  nature, 
the  direct  simulation  Monte  Carlo  method  works  with  far  fewer 
molecules  than  nature  does,  and  it  therefore  exhibits  consider¬ 
ably  greater  statistical  variation  in  its  macroscopic  predic¬ 
tions.  To  reduce  these  variations,  the  codes  are  run  repeatedly 
for  the  same  case,  increasing  the  statistical  base  from  which 
the  macroscopic  output  is  derived.  If  care  is  taken  to  use 
efficient  techniques,  such  as  described  in  this  report,  then 
useful  results  can  usually  be  obtained  with  a  modest  computa¬ 
tional  effort.  This  statement  must  be  tempered,  however,  by  a 
realization  of  the  convergence  rate  for  Monte  Carlo  sampling. 
Basically,  the  statistical  error  in  the  output  converges  as 
one  over  the  square  root  of  the  sample  size  (or  run  time) . 

Hence,  if  a  solution  looks  good,  but  the  user  decides  he  would 
like  one  more  significant  digit  (i.e.,  he  would  like  the  statis¬ 
tical  error  to  be  reduced  to  0.1  times  its  current  value)  it 
would  require  that  the  run  time  be  increased  by  a  factor  of  100! 
It  can  be  seen  that  the  desire  for  more  accuracy  can  very  quickly 
turn  the  most  efficient  code  into  a  money  gobbling  nightmare. 

When  using  a  Monte  Carlo  technique,  one  must  accept  some  statis¬ 
tical  scatter  in  the  output. 

8.2  Sampling  of  Instantaneous  Output  Quantities 

Instantaneous  output  quantities  are  those  which  are,  in 
principle,  derivable  from  an  instantaneous  "snapshot"  of  the 
solution.  These  quantities,  such  as  density,  temperature  and 
velocity,  can  be  determined  by  examining  the  molecular  state 
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vectors  at  a  particular  time  in  the  simulation.  The  code  pauses 
in  the  simulation  and  uses  the  molecular  state  vector  elements 
to  add  values  to  statistical  sums  appropriate  to  the  various 
cells  and  the  particular  time  that  it  passed.  It  then  proceeds 
with  the  simulation  until  the  next  sampling  time.  As  the  code 
goes  through  its  successive  runs,  it  stops  at  the  same  points  in 
the  simulation  every  time  and  adds  to  the  statistical  base  for 
the  sums.  The  items  listed  below,  with  their  statistical  defi¬ 
nitions,  are  selectable  as  output  requests  in  the  SSI  codes. 
Summations  are  made  over  all  applicable  simulated  molecules, 

which  includes  N  „  separate  runs. 
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Overall  Translational  Temperature 

T  -  _L_  (s  -  S3  +  j  *  S5 

3R0S1  V  2  S6 


where 


!i  -  Z  wi 


(101) 


(102) 


S2  =  52  Wimi  (Vli  +  V2i  +  v3i) 


(103) 


S3  =  X)  Wimivli 


(104) 


S4  *  X  WimiV2i 


(105) 


S5  *  H  Wiraiv3i 


(106) 


S6  H  Vi 


(107) 


Translational  Temperature  in  j'th  Direction 
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where 


S7  -  Z  Vivji 
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Internal  Mode  Temperature 
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with  the  exception  of  Eq.  (99) ,  all  of  the  above  quantities  can 
also  be  defined  and  calculated  for  any  specified  species.  The 
sums  are  the  same  except  that  only  molecules  of  that  species 
are  considered.  Before  printing  output  quantities,  they  are 
always  transformed  to  standard  dimensions  from  the  int^tnal 
scales. 


8.3  Sampling  of  Time  Averaged  Quantities 

Some  additional  quantities  of  interest  are  not  sampled  at 
a  separate  sampling  time  as  described  above,  but  rather  as  the 
simulation  evolves.  Examples  of  such  quantities  are  collision 
rates,  reaction  rates,  mean  velocities  between  molecules,  etc. 
In  general,  these  quantities  depend  on  the  relative  state  of 
more  than  one  type  of  molecule,  and  they  are  by  their  nature 
expressed  as  average  values  over  a  finite  time  interval.  The 
formulas  for  calculating  these  quantities  are  little  more  than 


event  counters,  and  will  not  be  included  here.  The  following 
quantities  are  currently  available  as  output: 

•  Mean  Relative  Velocity  Between  any  Two  Species 

•  R.M.S.  Deviation  of  Mean  Relative  Velocity 
Between  any  Two  Species 

•  Mean  Product  of  Cross  Section  Times  Relative 
Velocity  Between  any  Two  Species 

•  Collision  Rate  Between  any  Two  Species 

•  Reaction  Rate  for  any  Chemical  Reaction. 

The  sampling  for  all  of  these  quantities  occurs  in  the  colli¬ 
sion  simulation  routines.  As  pairs  are  considered  as  possible 
collision  partners,  statistics  are  kept,  if  necessary,  to  gener¬ 
ate  the  first  three  quantities.  Statistics  on  collisions  and 
reactions  are  kept  as  they  occur. 

9.  APPLICATION  TO  THE  FREEJET  EXPANSION  WITHIN  A 

MASS  SPECTROMETER 

9.1  Motivation  for  Studying  Major  Species  Freejet 

The  freejet  expansion  of  the  major  (neutral)  species  through 
a  sonic  orifice  into  a  near  vacuum  is  a  classic  problem  repre¬ 
senting  the  zeroth  order  solution  to  the  problem  at  hand.  The 
numerical  dominance  of  the  major  species  assures  that  their 
distribution  in  phase  space  (i.e.,  their  density  and  velocity 
distributions)  will  be  negligibly  affected  by  the  minor  species 
of  interest  (e.g.,  ion  clusters,  HjO,  etc.).  Hence,  for  a  given 
atmospheric  pressure  and  orifice  geometry  the  freejet  expansion 
of  the  major  species  can  be  calculated  once  and  for  all.  The 
physical  processes  of  interest  (ion  acceleration  due  to  electric 
fields,  agglomeration,  fragmentation,  etc.)  can  then  be  handled 
by  considering  the  minor  species  to  be  traveling  within  a  known 


phase  space  distribution  of  major  species.  Considerable  concep¬ 
tual  and  computational  simplifications  result  from  the  recogni¬ 
tion  that  the  free jet  expansion  is  uncoupled  from  the  ionic 
motion. 

9.2  Physical  Description  of  the  Free jet  Expansion 

The  major  features  of  a  freejet  created  by  expanding  air 
through  a  sonic  orifice  into  a  near  vacuum  are  illustrated  in 
Figure  1.  For  the  cases  of  interest  (20-40  km  altitude,  orifice 
diameter  =  0.04  cm)  the  orifice  is  always  considerably  larger 
than  an  ambient  mean  free  path,  so  the  initial  portion  of  the 
expansion  is  best  described  in  terms  of  continuum  fluid 
mechanics.  The  orifice  forms  a  sonic  throat  in  which  the  flow 
is  accelerated  to  a  Mach  number  of  unity.  If  the  orifice  has 
a  finite  thickness  ("t"  in  Figure  1),  then  a  boundary  layer 
forms  on  the  edge  of  the  orifice,  restricting  the  flow  somewhat 
from  that  which  would  be  predicted  via  one-dimensional  inviscid 
theory.  The  mass  flow  through  the  orifice  can  be  represented 


2  _ _  ,  ,  <<Y  +  1)/[2(y  -  1)]) 

m  =  CD  ^  YP0P0  [2/(y  +  U]  (112) 

where  rQ  is  the  orifice  radius,  y  is  the  ratio  of  specific 
heats  in  the  gas  and  PQ  and  pQ  are  the  stagnation  pressure  and 
density,  respectively.  The  discharge  coefficient,  CD,  is  a 
corrective  factor  to  bridge  the  gap  between  the  ideal  and  real 
worlds.  CD  is  less  than  unity  due  to  the  presence  of  the  ori¬ 
fice  boundary  layer  (viscosity  influence)  and  two  dimensional 
flow  effects. 

^Shapiro,  A.H.,  The  Dynamics  and  Thermodynamics  of  Compressible 
Fluid  Flow,  the  Ronald  Press  Co.,  New  Vork,  1953,  pp.  73-105. 


As  the  gas  passes  through  the  orifice,  an  expansion  fan  is 
formed  on  the  inner  edge  of  the  orifice  and  spreads  out  into  the 
flowfield.  The  expansion  fan  is  initially  describable  in  terms 
of  Prandtl-Meyer  theory,  but  is  altered  from  that  form  as  it 
proceeds  away  from  the  orifice  by  axisymmetric  (as  opposed  to 
planar)  flow  effects. 

The  initial  portion  of  the  supersonic  expansion  is  well 
suited  for  calculation  via  the  method  of  characteristics  (MOC),6 
although  the  interaction  of  the  boundary  layer  with  the  expan¬ 
sion  fan  is  difficult  to  treat  analytically.  As  the  flow  pro¬ 
ceeds  away  from  the  orifice  the  rotational  and  randomized 
translational  energy  is  transferred  to  directed  kinetic  energy 
of  the  flow.  The  flow  expansion,  and  consequent  reduction  in 
collision  frequency,  causes  the  gas  to  depart  from  thermal 
equilibrium.  A  residual  rotational  energy  remains  in  the  gas 
after  collisions  cease,  and  this  is  often  characterized  by  a 
rotational  temperature.  (The  populated  rotational  states  do 
not,  in  general,  adhere  to  a  Boltzmann  distribution;7  but  this 
is  of  no  great  consequence  to  the  present  investigation.)  Fur¬ 
thermore,  the  expansion  is  characterized  by  distinct  parallel 
and  perpendicular  temperatures  representing  residual  randomized 
kinetic  energy  in  directions  parallel  to  and  perpendicular  to 
the  mean  flow  direction. 

The  portion  of  the  jet  at  large  angles  from  the  symmetry 
axis  is  dominated  by  the  slower,  more  easily  turned,  gas  from 
the  orifice  boundary  layer.  Also,  the  lighter  molecular  weight 

7 

Sharafudinov,  R.G.  and  Skovorodko,  P.A.,  "Rotational  Level 
Population  Kinetics  in  Nitrogen  Freejets,"  Proceedings  of  the 
12th  International  Symposium  on  Rarefied  Gas  Dynamics,  AIAA 
Press,  1980. 


species  become  more  concentrated  at  large  angles  trom  the  sym¬ 
metry  axis  leaving  the  central  portion  of  the  jet  more  concen- 

o 

trated  in  species  with  heavier  molecular  weights. 

9.3  Potential  Calculational  Approaches 

For  the  most  part  the  departure  of  the  free jet  from  local 
thermodynamic  equilibrium  marks  the  end  of  the  region  of  valid¬ 
ity  for  continuum  techniques  such  as  the  method  of  character¬ 
istics.  (An  exception  to  this  is  a  recent  paper  by  Labowski, 

Ryali,  and  Fenn  using  the  so-called  nonequilibrium  method  of 

.  9 

characteristics  on  a  rotationally  relaxing  free^et.  However, 
the  method  does  not  apply  when  the  translational  modes  go  out 
of  equilibrium.  Since  the  region  between  rotational  and  trans¬ 
lational  nonequilibrium  is  usually  small,  the  method  is  inappro¬ 
priate  in  the  present  study. ) 

The  only  presently  available  computational  technique  which 
is  capable  of  describing  the  critical  molecular  processes 
occurring  in  the  latter  portions  of  the  expansion  is  the  direct 
simulation  Monte  Carlo  method.  There  is  no  choice  to  be  made 
for  the  calculational  technique  to  be  used  in  the  nonequilibrium 
portion  of  the  expansion;  the  only  question  is  whether  the  Monte 
Carlo  method  should  be  used  for  the  entire  solution  region,  or 
whether  it  should  be  wedded  to  a  MOC  calculation  of  the  equil¬ 
ibrium  portion.  The  answer  to  this  question  depends  largely  on 
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Bird,  G.A.,  "Breakdown  of  Continuum  Flow  in  Freejets  and  Rocket 
Plumes,"  Proceedings  of  the  12th  International  Symposium  on 
Rarefied  Gas  Dynamics,  AIAA  Press,  1980. 
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Labowski,  M.  Ryali,  S.,  and  Fenn,  J.B.,  "Flowfield  Calculations 
in  Nonequilibrium  Freejets  by  the  Method  of  Characteristics," 
Proceedings  of  the  12th  International  Symposium  on  Rarefied 
Gas  Dynamics,  AIAA  Press,  1980. 


how  large  the  equilibrium  region  is  expected  to  be.  (It  should 
be  stressed  that  the  Monte  Carlo  flowfield  is  valid  for  equil¬ 
ibrium  continuum  flow  —  but  it  is  less  efficient  computationally.) 


The  breakdown  of  translational  and  rotational  equilibrium 
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has  been  considered  in  some  detail  by  Bird,  '  and  it  has  been 
observed  to  correlate  well  with  the  paparameter  P  given  by 


1  D (&np) 
V  Dt 


(113) 


where  p  is  the  fluid  density,  v  is  the  collision  frequency  and 
D/Dt  represents  the  substantial  derivative  with  respect  to  time. 
The  breakdown  of  equilibrium  has  been  observed  to  correlate 
quite  well  with  a  P  value  of  approximately  0.05. 

Figure  2  shows  contours  of  constant  P/Kn  as  calculated  in 
Ref.  8  for  the  expansion  of  a  diatomic  gas  from  an  orifice  into 
a  vacuum.  In  this  figure  the  Knudsen  number,  Kn,  is  defined  as 
the  ratio  of  the  mean  free  path  at  the  orifice  to  the  orifice 
radius.  The  calculation  was  performed  via  the  method  of  charac¬ 
teristics,  so  the  orifice  Mach  number  was  taken  to  be  1.1  to 
allow  the  calculational  grid  to  step  away  from  the  orifice. 

For  the  parameter  range  of  interest  (altitudes  of  20-40  km, 
orifice  radius  =  0.02  cm)  Figure  2  can  be  used  to  estimate  the 
position  at  which  the  onset  of  nonequilibrium  occurs,  and, 
hence,  the  maximum  range  of  validity  for  a  MOC  calculation 
before  it  must  be  joined  to  a  Monte  Carlo  calculation.  The 
conclusion  is  that  a  MOC  calculation  cannot  be  used  beyond  two 
orifice  radii  from  the  orifice  at  20  km  and  essentially  cannot 
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Figure  2.  Contours  of  Mach  Number  and  Continuum 
Breakdown  Criterion  in  the  Free jet  Expansion  of  a 
Diatomic  Gas.  (From  Reference  8.) 


be  used  at  all  at  40  km  altitude.  This  very  limited  range  of 
validity  for  a  MOC  calculation  strongly  suggests  that  the  pre¬ 
ferred  calculational  approach  is  simply  to  use  a  Monte  Carlo 
method  for  the  entire  flowfield. 


9.4  Cell  Mesh  Selection 

The  selection  of  grid  geometries  for  fluid  mechanic  calcu 
lations  can  generally  be  regarded  as  more  of  an  art  than  a 
science.  Considerations  in  the  selection  of  a  grid  are: 

•  The  grid  should  be  as  simple  as  possible,  since 
the  program  must  repeatedly  decide  what  cell  mole¬ 
cules  are  in  as  they  move.  If  this  required  the 
solution  of  a  complex  equation,  the  entire  pro¬ 
gram  would  run  significantly  slower  than  if  the 
cell  can  be  determined  easily. 


•  The  grid  should  concentrate  cells  where  gradients 
are  the  largest,  so  that  the  least  number  of  total 
cells  (and  molecules)  are  needed  to  obtain  an 
accurate  solution. 

•  The  grid  should  provide  output  where  it  is  required, 
with  the  resolution  that  is  desired  for  the  answer 
of  interest. 


The  cell  structure  which  was  introduced  in  the  first 
quarterly  technical  report  was  based  mainly  on  the  second  con¬ 
sideration  listed  above.  It  concentrated  cells  in  the  vicinity 
of  the  Prandtl -Meyer  discontinuity,  where  the  flow  gradients 
are  largest.  It  suffers  disadvantages,  however,  in  that  spur¬ 
ious  small  cells  occur  on  the  axis  and  cell  determination  for 
a  molecule  is  somewhat  difficult.  As  a  first  attempt,  therefore, 
a  simpler  cell  structure  was  devised.  The  cell  structure  is 
illustrated  schematically  in  Figure  3,  and  is  characterized  by 
the  following  relations: 


•  The  orifice  radius  is  divided  up  evenly  into  a 
specified  number  of  grid  spacings.  This  number 
can  be  varied  to  achieve  convergence. 

•  Above  the  orifice,  at  the  first  axial  location, 
there  are  some  additional  cells  of  the  same  size 
to  help  describe  the  expansion  that  occurs  around 
the  orifice  edge.  The  number  of  cells  above  the 
orifice  can  also  be  varied  to  achieve  convergence 
along  the  flow  centerline  (the  region  of  import¬ 
ance  for  the  present  investigation) . 

•  The  cells  are  constructed  between  a  series  of 
planes  which  are  perpendicular  to  the  jet  axis. 

The  spacing  between  successive  planes  increases 
as  they  proceed  further  from  the  orifice,  since 
the  flow  gradients  will  decrease  in  this  direction. 
The  rule  for  the  increase  of  the  spacing  between 
the  planes  is  given  by  the  equation: 


ro  +  zj 


A 


constant 


9 


(114) 


RADIAL  DIRECTION 


I 


s 


i 


['• 


where  Zj  is  the  location  of  the  j'th  plane  and 
DZj  is  the  spacing  between  it  and  the  previous 
plane  (i.e.,  DZj  =  Zj  -  Zj_j),  and  rQ  is  the 
orifice  radius.  This  relation  can  be  solved  to 
give  the  location  of  Zj  directly  as: 


This  relation  is  easily  inverted  to  find  the 
axial  cell  location  of  a  particular  molecule. 

•  The  cells  are  selected  to  have  a  radial  increment 
equal  to  the  axial  increment,  and  the  number  of 
cells  per  axial  location  is  kept  constant.  Hence, 
as  the  cells  grow  larger  in  the  axial  direction, 
they  automatically  keep  pace  in  the  radial  direction. 


9.5  Sample  Calculational  Results 

Initial  calculations  were  performed  for  atmospheric  condi¬ 
tions  appropriate  to  an  altitude  to  30  km  (ambient  number 
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density  of  3.83  x  10  molecules/cm  ,  ambient  temperature  of 
227  K)  for  a  mixture  of  79%  nitrogen  and  21%  oxygen.  The  flow 
was  expanded  isentropically  to  a  Mach  number  of  unity  at  the 
orifice,  and  the  flow  profile  was  assumed  uniform  across  the 
orifice.  As  was  discussed  previously,  this  is  recognized  to 
represent  an  idealization;  but  it  seemed  like  a  reasonable 
place  to  start  the  calculations. 

The  first  question  that  was  addressed  was  to  determine  what 
type  of  grid  was  required  to  achieve  convergence  of  the  solution. 
Figure  4  shows  plots  of  axial  number  density  obtained  from  two 
separate  grid  spacings,  for  otherwise  identical  runs.  The 
coarser  grid  divided  the  orifice  radius  into  four  separate  cells 
(which  means,  of  course,  that  the  diameter  was  divided  into 
eight  cells)  and  had  two  additional  cells  above  the  orifice. 

The  finer  grid  had  twice  as  many  cells  within  and  above  the 
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orifice  radius.  The  calculations  were  carried  out  for  approxi¬ 
mately  3.5  orifice  diameters  downstream  (0.15  cm),  by  which 
distance  the  flow  was  substantially  supersonic.  This  is  impor¬ 
tant  since  it  means  that  the  neglect  of  molecules  traveling 
upstream  into  the  solution  region  is  quite  well  justified. 

9.6  Discussion  of  Free jet  Calculations 

The  first  conclusion  to  draw  from  these  preliminary  calcu¬ 
lations  is  that  the  program  is  working,  and  the  results  look 
reasonable.  It  can  be  seen,  however,  that  the  coarse  grid 
spacing  is  not  adequate  to  achieve  a  convergence  of  the  solu¬ 
tion,  and,  on  the  basis  of  these  results,  it  cannot  be  said 
whether  or  not  the  fine  grid  spacing  is  converged.  Further 
calculations  will  clarify  this  issue  and  others  regarding  mesh 
convergence.  It  is  also  desirable  to  know  how  many  cells  are 
required  above  the  orifice  to  achieve  accurate  axial  profiles. 

If  the  required  number  of  cells  to  achieve  mesh  convergence 
for  the  cell  structure  outlined  in  this  section  is  too  large, 
then  it  may  well  be  necessary  to  adopt  another  structure  such 
as  the  one  outlined  previously.  It  is  expected  that  the  ideal 
orifice  considered  here  represents  a  worst  case,  since  when  an 
actual  flow  profile  is  introduced  across  the  orifice,  the  ori¬ 
fice  boundary  layer  will  lessen  the  importance  of  the  Prandtl- 
Meyer  discontinuity,  and  therefore  the  importance  of  small  cells 
in  that  region. 

Rotational  relaxation  was  calculated  for  these  test  cases, 
and  a  significantly  higher  rotational  than  translational  temper¬ 
ature  (106  K  versus  25  K)  was  predicted.  These  numbers  have  no 
quantitative  significance,  since  no  attempt  was  made  to  input  a 
proper  rotational  relaxation  collision  number.  However,  it  is 
noteworthy  that  the  code  is  capable  of  predicting  this  physically 
real  behavior. 


9.7  Real  Orifice  Flow  Effects 
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In  principle,  it  would  be  possible  to  compute  the  flow  from 
the  undisturbed  atmosphere  into  the  mass  spectrometer  in  a  single 
solution.  In  such  a  solution,  the  boundary  conditions  applied 
would  correspond  to: 

1)  An  undisturbed  atmosphere  bounding  the  appropriate 
portion  of  the  solution  region. 

2)  A  wall  interaction  representing  the  molecular 
collisions  with  the  mass  spectrometer  casing  as 
well  as  the  sides  of  the  orifice. 

3)  Possible  downstream  conditions  if  the  solution  is 
carried  to  the  point  where  interaction  with  the 
background  gas  within  the  mass  spectrometer  is 
important . 

Although  wall  boundary  conditions  do  have  a  degree  of 
approximation  that  reflects  our  incomplete  understanding  of  the 
gas-surface  interaction  process,  the  above  boundary  conditions 
are  relatively  well  characterized.  However,  the  above  procedure 
does  involve  a  substantial  amount  of  computational  effort  aimed 
at  the  flow  outside  of  the  mass  spectrometer,  since  that  flow  is 
usually  characterized  by  a  relatively  small  mean  free  path. 

Since  the  primary  interest  of  the  present  study  is  to  describe 
the  flow  within  the  mass  spectrometer,  it  is  natural  to  try  to 
reduce  the  size  of  the  solution  region  while  retaining  the 
essential  portion  of  the  flow.  The  obvious  way  to  do  this  is  to 
approximate  the  flow  at  the  orifice,  and  merely  compute  the  flow 
within  the  instrument. 

There  are  several  reasons  to  believe  that  such  a  procedure 
would  be  fruitful.  If  the  orifice  were  a  smooth  walled  converging- 
diverging  nozzle,  then  the  flow  would  expand  to  a  Mach  number  of 
unity  at  the  orifice  plane,  and  it  would  be  known  exactly  (except 
for  boundary  layer  and  two  dimensional  effects)  irrespective  of 
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the  flow  downstream  of  the  orifice.  Although  the  actual  orifice 
has  a  sharp  edge#  intuition  implies  that  the  downstream  flow  can 
have  at  best  a  small  effect  on  the  orifice  flow,  so  it  would 
seem  plausible  to  specify  the  orifice  flow  a  priori.  Further¬ 
more,  the  present  study  is  concerned  mainly  with  the  flow  along 
the  centerline  of  the  resulting  jet,  where  boundary  layer  and 
two  dimensional  effects  should  be  minimized. 

Although  it  would  seem  to  be  a  classic  problem,  the  expan¬ 
sion  from  a  uniform  state  through  a  sharp  edged  orifice  into  a 

* 

vacuum  has  apparently  not  been  the  object  of  intense  study. 
Liepmann^1  considered  the  problem  in  sonie  detail,  and  indicated 
that  the  flow  in  the  plane  of  the  orifice  is  subsonic  along  the 
centerline,  going  supersonic  and  then  back  to  sonic  as  the  edge 
of  the  orifice  is  approached.  Little  else  was  said  in  this  or 
any  other  available  reference  about  the  actual  flow  in  the  ori¬ 
fice  plane.  Measurements  were  presented,  but  they  were  limited 
to  ^ne  discharge  coefficient  (the  ratio  of  ideal  to  actual  mass 

flow),  which  was  subsequently  measured  by  Smetana,  Sherrill  and 
12 

Schort  to  greater  accuracy.  The  results  of  the  latter  study 
are  presented  in  Figure  5,  which  has  been  replotted  to  obtain 
a  more  convenient  form.  (As  presented  in  Ref.  12  both  axes 
involved  the  unknown  mass  flow,  so  an  iteration  was  required  to 
determine  the  discharge  coefficient.  This  problem  is  removed 
when  the  Reynold’s  number  is  defined  in  terms  of  available 
reference  quantities  rather  than  the  unknown  mass  flow.  The 
two  Reynold's  numbers  differ  only  by  a  factor  of  the  discharge 
coefficient. ) 

* ^Liepmann ,  H.W. ,  "Gaskinetics  and  Gasdynamics  of  Orifice  Flow," 
Journal  of  Fluid  Mechanics,  Vol.  10,  No.  5,  1961,  pp.  65-79. 

12 

Smetana,  F.O.,  Sherrill,  W.A.,  and  Schort,  D.R.,  "Measurements 
of  the  Discharge  Characteristics  of  Sharp-edged  and  Round-edged 
Orifices  in  the  Transition  Regime,"  Proceedings  of  the  5’th 
Rarefied  Gas  Dynamics  Symposium,  Vol.  2,  Academic  Press,  1967. 
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Figure  5.  The  Discharge  Coefficient  for  a  Sharp 
Edged  Orifice  as  Measured  by  Smetana, 
Sherrill  and  Schort.O)  it  is  plotted 
here  as  a  function  of  Reynold's  Number 
based  on  sonic  conditions  and  the  orifice 
diameter. 


Use  of  the  discharge  coefficient  presented  in  Figure  5 
accounts  for  the  major  effect  of  a  sharp  edged  orifice  as  opposed 
to  a  smooth  nozzle,  and  it  is  apparently  the  only  available 
effect  for  which  measurements  exist.  The  initial  approach  will 
be  to  adjust  the  orifice  mass  flow  to  reflect  the  results  of 
Figure  5,  but  the  mean  flow  velocity  will  still  be  assumed  to 
start  at  sonic  conditions.  (Individual  molecules,  of  course, 
will  have  a  distribution  of  velocities.  See  Section  7  for  a 
complete  description  of  the  handling  of  boundary  conditions.) 

Any  error  incurred  via  this  procedure  should  have  little  effect 
on  the  centerline  solution  of  interest. 

9.8  Handling  of  Minor  Species 

As  discussed  in  Section  6,  a  procedure  has  been  devised 
at  SSI  to  handle  sampling  of  collisions  simultaneously  between 
an  arbitrary  number  of  species.  For  the  mass  spectrometer  case, 
however,  it  may  be  advantageous  to  recast  this  procedure  in 
terms  of  major  and  minor  species.  The  solution  for  the  major 
species  can  be  carried  out  once  and  for  all  for  a  given  altitude 
and  instrument  design.  Then,  since  the  minor  species  are  too 
numerically  insignificant  to  affect  the  major  species  flow,  it 
is  possible  to  consider  several  possible  atmospheric  concentra¬ 
tions  of  various  minor  species  in  order  to  obtain  a  best  fit 
with  the  data.  In  order  for  this  procedure  to  work,  of  course, 
there  has  to  be  a  method  to  go  back  and  calculate  minor  species 
distributions  given  the  major  species  solution.  It  is  the  pur¬ 
pose  of  this  section  to  present  that  method. 

In  such  a  minor  species  solution,  it  is  only  necessary  to 
define  a  density,  velocity  and  temperature  in  each  cell  for  each 
of  the  major  species.  (Actually  a  different  effective  tempera¬ 
ture  would  be  defined  for  each  direction.)  Then,  for  purposes 
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of  sampling  collisions  between  major  and  minor  species,  major 
s;pecies  molecules  can  be  generated  from  the  proper  distribution 
as  needed.  These  molecules  need  not  be  permanently  stored,  nor 
is  it  necessary  to  advance  them  along  their  trajectories. 

The  major  computational  advantage,  however,  is  the  removal 
of  a  necessity  for  calculating  collisions  of  major  species  with 
themselves.  Since  collision  sampling  is  the  most  time  consuming 
element  of  a  Monte  Carlo  simulation,  and  the  collisions  between 
major  species  are  the  most  prevalent,  there  is  room  for  a  sub¬ 
stantial  increase  in  computational  speed.  Within  the  framework 
described  in  this  report,  this  can  be  accomplished  by  defining 
separate  global  time  counters  for  major-minor  and  minor-minor 
collisions.  Each  global  time  counter  is  defined  as  a  weighted 
average  of  the  time  counters  for  each  collision  class  that 
falls  within  its  domain.  The  weighting  factors  are  selected 
to  achieve  tractable  forms  for  the  two  time  counters. 

9.8.1  Time  Counter  for  Collisions  Between  Two  Minor 
Species 

Ordinarily,  collisions  between  two  minor  species  would  be 
considered  unimportant  in  comparison  to  collisions  between  major 
and  minor  species.  However,  for  the  present  study  there  are 
processes  such  as  agglomeration  which  may  not  be  negligible 
(due  to  very  large  cross  sections)  even  though  they  involve 
two  minor  species.  The  handling  of  minor  species  self-collisions 
can  be  achieved  exactly  as  described  in  Section  6,  with  the 
adjustment  that  total  number  density  is  taken  to  refer  to  total 
minor  species  number  density.  The  time  counter  increment  is 
then  just  as  given  in  Eq.  (76). 


9.8.2  Time  Counter  for  Collisions  Between  Major  and 
Minor  Species 

A  separate  global  time  counter  is  easily  defined  for  major- 
minor  collisions.  If  P  denotes  the  number  of  major  species  and 
p  denotes  the  number  of  minor  species,  then  the  major-minor 
collision  time  counter,  tg,  can  be  defined  via 
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which  is  the  obvious  extendion  of  Eq.  (72) . 
choice  for  ^  proves  to  be 
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The  convenient 
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>In  the  above  relations,  an  i  subscript  is  now  taken  to  refer 
to  a  major  species,  and  a  j  subscript  is  taken  to  refer  to  a 
minor  species.  In  Section  6  the  subscripts  make  no  such  dis¬ 
tinction.)  Substituting  Eq.  (117)  into  Eq.  (116),  and  carrying 
out  the  summation  in  the  denominator  yields  an  expression  for 
the  major-minor  global  time  counter  in  terms  of  the  time 
counters  for  the  individual  collision  classes.  The  appropriate 
increment  in  the  global  time  counter  to  make  when  simulating  a 
collision  is 
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where  n^  and  nj  refer  to  the  total  number  densities  of  all  major 
and  minor  species,  respectively,  V  is  the  cell  volume,  and  o^j 
is  the  mutual  cross  section  for  the  species  colliding  at  a 
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relative  velocity  cr>  The  time  counter  increment  represented 
by  Eq.  (118)  is  much  larger  than  the  increment  would  be  for 
the  same  collision  (Eq.  (76))  if  minor  species  were  not  being 
handled  separately.  This  is  the  quantitative  realization  of 
the  general  principle  that  collision  sampling  goes  a  lot  more 
quickly  if  major-major  collisions  need  not  be  considered. 
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